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FOREWORD 


This  training  guide  explains  the  fundamentals  and  foundations  of  the  long-range  ballistic 
missile  (LRBM)  weapon  systems  launched  from  moving  bases  such  as  submarine  or  aircraft. 
Since  the  land-based  missiles  are  essentially  no  different  from  the  missiles  launched  from  the 
moving  bases  except  that  the  bases  are  fixed  on  the  earth  (without  rotational  and  linear  motions), 
the  theory  and  analyses  explained  in  this  training  guide  would  be  equally  applicable  to  all  three 
classes  of  missiles — land-based,  sea/undersea-based,  or  airborne  missiles,  which  use  inertial 
navigation  system. 

It  is  written  mainly  for  the  training  of  new  technical  members  using  only  introductory 
physics  and  calculus.  Nonetheless,  it  manages  to  maintain  the  robustness  of  theory  and  the  rigor 
of  mathematical  derivation.  For  this  reason,  it  is  recommended  to  veteran  technical  members  as 
well,  who  might  find  a  segment  or  two  that  would  shed  a  new  light  on  old  subjects. 
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1.0  INTRODUCTION 


This  publication  is  a  training  guide  explaining  fundamentals  and  foundations  of  the  long- 
range  ballistic  missiles  (LRBM)  from  the  viewpoint  of  the  navigation  and  guidance.  The  missiles 
may  be  launched  from  either  fixed-base  on  land  or  from  moving-base  such  as  aircraft  or 
submarine.  The  articles  are  written  at  such  a  level  that  anyone  with  freshman  physics  and  calculus 
would  not  have  much  difficulty  in  understanding  them.  Readers  of  this  publication  are  cautioned 
that  the  analysis  of  the  actual,  real-world  systems  are  much  more  complicated  than  those  described 
here.  As  the  title  implies,  these  articles  are  merely  "fundamentals  and  foundations." 

Section  2.0  lists  the  symbols  and  notations  used  in  this  report.  It  includes  an  operator  P  ( • ) 
for  -j-  ( • ) .  This  is  handy  because  the  differentiation  with  respect  to  time  in  the  earth-fixed  E-frame 

may  simply  be  expressed  as  PE  ( •) .  Section  3.0  describes  the  frames  of  reference  relevant  to  this 

report.  Section  4.0  explains  implementation  of  Newton's  second  law  of  motion  for  the  purpose  of 
inertial  navigation.  Since  all  LRBMs  use  the  center  of  the  earth  as  the  origin  of  the  inertial 
reference  frame.  Section  5.0  deals  with  the  errors  incurred  in  the  earth-centered  inertial  frame. 
Section  6.0  explains  that  the  effect  of  angular  rotations  are  dependent  on  the  sequence  in  which 
they  are  executed,  unlike  vector  additions.  Theory  of  inertial  navigation  involves  rotation  of  one 
reference  frame  relative  to  another.  The  theorem  of  Coriolis,  which  is  discussed  in  Section  7.0,  is 
essential  to  understand  the  relationship  between  the  two  rotating  frames  relative  to  each  other.  The 
matrix  form  of  the  theorem  of  Coriolis  is  discussed  in  Section  8.0.  Section  9.0  discusses  the 
rotational  form  of  Newton's  second  law  of  motion,  which  is  essential  for  the  understanding  of  the 
operation  of  gyros,  gyro-accelerometers,  and  stable  platforms  of  guidance  system  of  LRBMs. 
Section  10.0  describes  the  prototype  linear  accelerometer,  which  serves  as  a  tutorial  for 
understanding  the  pulse  integrated  pendulous  accelerometers  (PIPA),  and  pendulous  integrating 
gyro  accelerometers  (PIGAs),  is  described  in  Section  11.0.  Section  12.0  discusses  the  rotational 
motion  of  the  missile  body  based  on  some  simple  assumptions,  which  allows  us  to  avoid 
complicated  mathematics. 

Based  on  theoretical  and  mathematical  tools  developed  in  previous  sections,  the  theory  of 
practical  gyros  are  explained  in  Section  13.0,  the  single  degree-of-freedom  (SDOF)  gyro  in  Section 
14.0,  the  stable  platform  of  inertial  measurement  unit  in  Section  15.0,  and  the  PIGA  in  Section 
16.0. 


Moving  from  components  to  the  system  level,  the  derivation  and  implementation  of  the 
equation  of  motion  for  the  inertial  navigation  is  described  in  Section  17.0.  The  determination  of 
the  initial  position  and  initial  velocity  to  be  used  in  the  integrations  by  the  guidance  computer  of  the 
Newton's  second  law  of  motion  is  described  in  Section  18.0.  Moving  to  the  rocket  steering 
scheme  during  powered  stages,  the  so-called  cross  product  steering  is  described  in  Section  19.0. 

Toward  the  end  of  the  powered  flight,  the  guidance  system  of  the  missile  tries  to  correct  its 
error,  in  a  statistical  sense,  by  viewing  a  predesignated  star.  The  weighing  matrix,  also  called  the 
W  matrix,  which  is  used  in  the  statistical  computation  for  this  purpose,  is  derived  in  Section  20.0. 
The  topic  of  Schuler  tuning  is  discussed  concisely  in  Section  21.0.  Schuler  tuning  is  conceptually 
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important  because  it  deals  with  the  prevention  of  platform  misorientation  caused  by  the  acceleration 
in  the  earth's  gravitational  field.  Finally,  the  correlated  velocity,  which  is  the  required  velocity  of 
the  reentry  vehicle  to  hit  the  target  on  the  rotating  earth  by  free  falling  with  specified  time-of-flight 
(TOF),  is  discussed  in  Section  22.0. 


1-2 


NSWCDD/MP-95/87 


2.0  LIST  OF  SYMBOLS  AND  NOTATIONS 

Rpq  Displacement  vector  from  point  P  to  point  Q.  Similarly  Vpq  is  the  velocity  vector  of  Q 

relative  to  P,  etc. 

Wab  Angular  velocity  of  Frame  B  relative  to  Frame  A. 

VA  Any  vector,  (e.g.,  R  or  W,  etc.)  with  components  expressed  in  Frame  A. 

P()=> 

£(•)■/(•>* 

dA 

PA  ( •) s  — j^-  ( •)  Time  derivative  with  respect  to  (WRT )  Frame  A,  (i.e.,  the  increment  observed  in 
Frame  A). 

Thus, 


y 


p»R  =  p«(p«R)  =  p(pRa)  = 

PaPbR  =  Pa(PbR) 


dt2 


The  advantages  of  these  notations  will  become  obvious  in  the  second-half  of  this  report. 
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R°=  y  andW°=wy  , 


then  the  matrix  equivalent  of  the  vector  W°  ,  denoted  by  is  given  by 


'  0  -w,  wy 

wdk  w  0  -  w 

T  AB 

rwy  wx  0 


so  that  [W°J  ]r«W^b  xRd  or 


0  -wz  wY  Yx^  i  j  k  [  wyz- wzy 

wz  0  -wx  I  y  <=>  wx  wY  wz  <=>  WzX-WxZ 

-wY  wx  0  Xz)  x  y  z  D  |wxy-wyX> 


C“  coordinate  transformation  matrix  from  Frame  A  to  Frame  B. 


Note: 


pNpB  _  pN 


CtCB=CA=I  identity  matrix. 
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3.0  FRAMES  OF  REFERENCE 


3.1  INERTIAL  (I)-FRAME 

The  frame  in  which  Newton’s  laws  of  motion  (specifically  the  second  law)  is  valid  is  called 
inertial  (I)-frame  by  definition.  In  all,  LRBM,  an  earth-centered  nonrotating  frame  in  the  inertial 
space,  is  accepted  as  a  (quasi)  I-frame.  Obviously,  since  the  center  of  the  earth  accelerates  around 
the  sun  in  an  elliptic  orbit,  the  earth-centered  I-frame  (denoted  by  IE)  is  not  an  1-frame  in  a  strict 
sense.  However,  the  error  incurred  in  this  practice  is  negligible  (this  is  shown  in  detail  under  a 
separate  heading  in  the  sequel). 

In  the  standard  IE -frame,  the  z-axis  is  collinear  with  the  earth’s  polar  (spin)  axis.  The  x- 
axis  and  y-axis  are  on  the  earth’s  equatorial  plane  (but  not  rotating  with  the  earth)  in  such  a  way  as 
to  form  an  orthogonal  triad  by  the  right-hand  rule.  This  frame  is  sometimes  referred  to  as  the 
"polar  earth-centered  inertial  frame." 


3.2  EARTH  (E)-CENTERED,  EARTH-FIXED  FRAME 

This  frame,  being  earth  fixed,  rotates  with  the  earth  relative  to  the  inertial  space.  The  x,  y, 
and  z  axes  are  defined  similar  to  the  standard  IE  frame,  except  x  and  y  axes  are  fixed  to  the  earth 
and  therefore  rotate  with  the  earth. 


3.3  NAVIGATIONAL  (N),  LOCAL-LEVEL  FRAME 

The  origin  of  this  frame  is  located  at  the  center  of  the  onboard  inertial  navigation  system 
(INS)  of  the  ship,  submarine,  or  aircraft  cruising  on  or  near  the  surface  of  the  earth.  The  x-axis 
points  to  the  north  and  y-axis  to  the  east,  both  on  the  local-level  plane,  which  is  tangent  to  the 
reference  ellipsoid.  The  z-axis  points  downward  perpendicular  to  the  tangential  plane.  The  z-axis 
may  or  may  not  coincide  with  the  local  vertical  (plumb  bob)  lim  because  the  earth’s  gravity  field  is 
not  uniform.  This  causes  the  deflection  of  the  vertical. 


3.4  GUIDANCE  (G)- SYSTEM  INERTIAL  MEASUREMENT  UNIT  (IMU)  FRAME 

The  origin  of  the  G-frame  is  located  at  the  center  of  the  missile  guidance  system’s  IMU, 
with  its  axes  parallel  to  the  earth-centered  1-frame.  Obviously,  since  the  missile  is  subject  to 
acceleration,  the  G-frame  cannot  be  an  I-frame.  However,  since  its  axes  are  parallel  to  the  axes  of 

the  I-frame,  the  coordinate  transformation  matrix  from  the  G-frame  to  the  1-ffame  denoted  by  C:G 
must  be  an  identity  matrix  or  Cq= I,  where  I  denotes  the  identity  matrix. 
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The  onboard  accelerometer  measures  the  acceleration  of  the  instrument  relative  to  the  origin 
of  the  I-frame,  not  relative  to  the  origin  of  the  G-frame  to  which  the  instrument  has  a  fixed  distance 
and  therefore  the  relative  acceleration  should  be  zero.  However,  since  the  axes  of  G-frame  is 
parallel  to  the  axes  of  1-frame,  the  accelerometer  output  (sensed  acceleration  and  the  velocity-gain 
caused  by  it  during  the  computation  interval)  should  have  the  same  values  whether  they  are 
expressed  in  I-frame  or  G-frame.  That  is, 

AV1  =  CqAVg  =  IAVC  =  AV° 


where 

AV°  =  velocity  gain  (AV)  determined  by  the  guidance  system  with  components  expressed  in 
G-frame,  and 

AV1  =  AV  with  components  expressed  in  I-frame. 

Thus 

ARI=AVIAt  =  AV°At 
where 

AR1  =  position  displacement  with  components  expressed  in  the  I-frame 

We  see  that  as  far  as  the  determination  of  the  missile’s  incremental  position  and  velocity  are 
concerned  the  G-frame  behaves  as  if  it  were  an  I-frame,  although  in  reality  it  is  not.  To  determine 
the  current  position  or  velocity  relative  to  the  I-frame,  the  incremental  position  or  velocity  obtained 
in  the  guidance  frame  is  added  by  the  guidance  computer  to  the  previous  position  or  velocity, 
referenced  to  the  origin  of  the  I-frame,  which  is  conventionally  the  center  of  the  earth  in  all  LRBM 
applications.  For  this  reason,  some  incorrectly  call  the  G-frame  “inertial,”  which  should  be 
avoided. 
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4.0  NEWTON’S  SECOND  LAW 
(FOR  THE  IMPLEMENTATION  OF  INS) 


According  to  Newton’s  second  law 

mPX  =  ZF  (4-1) 

Since  the  sum  of  externally  applied  forces  XF  has  to  be  the  sum  of  the  gravitational  force 
Fg  and  the  nongravitational  force  Fnq,  we  may  write  (1)  as  follows: 


pi2rip  =  — +  — =  g  +  f 
m  m 


(4-2) 


where  g  is  the  gravitational  force  per  unit  mass  and  f  is  the  nongravitational  force  per  unit  mass. 

By  convention,  f  is  called  simply  “specific  force,”  meaning  specific  (per  unit  mass) 
nongravitational  force. 


The  reason  we  separate  f=  from  g=-^  is  by  necessity,  not  by  choice,  because  the 

on-board  inertial  instrument  (accelerometer)  can  measure  only  f,  not  g.  This  is  explained  in  the 
section  of  the  accelerometer  later.  The  g  is  determined  by  the  onboard  computer,  based  on  the 
mathematical  model  (inverse  square  law)  as  a  function  of  the  position  from  the  center  of  the 
earth,  which  is  the  origin  of  the  I-frame  by  convention. 

Thus,  in  the  INS,  Newton’s  second  law  (4-2)  is  implemented  as  shown  in  Figure  4-1. 


00 

Y 

UJ 


4-1 


NSWCDD/MP-95/87 


5.0  ERRORS  IN  THE  EARTH- CENTERED  INERTIAL  FRAME  (IE  FRAME) 


For  the  purpose  of  estimating  the  magnitude  of  the  errors  incurred  by  treating  the  earth- 
centered  inertially  nonrotating  frame  (IE  frame)  as  the  I-frame,  consider  a  system  consisting  of  the 
sun  (S),  the  earth  (E),  the  moon  (M),  and  a  missile  (P)  near  the  earth  with  the  center  of  the  sun 
considered  as  the  origin  of  an  I-frame  as  shown  in  Figure  5-1  (Figure  not  to  scale;  the  rest  of  the 
universe  is  ignored). 


FIGURE  5-1.  THE  SUN,  MOON,  AND  EARTH  SYSTEMS 


The  forces  acting  on  missile  P  with  mass  m  are  the  gravitational  forces  by  the  sun,  the 
earth,  the  moon,  and  the  nongravitational  forces  of  thrust  and  aerodynamic  force  (if  in  the 
atmosphere). 

By  applying  Newton’s  second  law  along  with  the  law  of  gravitation  to  the  missile  at  P. 

'Ump  ^ng  =  Ksp 


GMm  GMFmTT  GM_mTT  _  _2, 

_ 5 T T _ L  T T  m  TT  i  E  

—  o  JS.SP  —  U.  PP 


idEP 


^SP 


VEP 


RL 

mp 


(5-1) 


m 

as 

mass  of  the  missile  at  P 

Ms 

= 

mass  of  the  sun  at  S 

Me 

= 

mass  of  the  earth  at  E 

Mm 

as 

mass  of  the  moon  at  M 

IIsp 

= 

unit  vector  along  SP 

Uep 

= 

unit  vector  along  EP 

Ump 

= 

unit  vector  along  MP 

G 

= 

gravitational  constant 
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Since  R^p  —  Rse  +  Rep 


mP2Rsp  =  mP,2  Rse  +  mP2REP 
Dividing  by  (1),  by  m  and  using  (2): 

GMStt  GMEtt  .  F, 


•uSP-^uFP- 


p2  ^-SP  p2  iiEP  p2 
R-Sp  R-BP  Rmp 


Ump  +  ^  =  Pi2Rse+Pi2Rep 
m 


Now,  applying  Newton’s  second  law  to  the  earth: 

Til  p*D  _  GMsMe  GMmMe 

ME*IKSB  p2  —SB  p  2  —ME 

^SB  ^ME 


(5-2) 


(5-3) 


(5-4) 


or 


2R-  =  - 


PjR 


SE 


GM»n 


,2  — SE 

‘‘SE 


ME 


‘‘ME 


(5-5) 


Substituting  (5)  into  (3)  for  P^R 


I  SE 


2r»  GM 

Pi  Rep  p2 

kep 


^Ep+^+Ags+AgM 


m 


(5-6) 


where 


Ags  =  -GMS 


1 


1 


^2-USp  2  USE 

VKSP  ^SE  J 


(5-7) 


is  the  difference  between  gravitational  forces  per  unit  mass  by  the  sun  on  the  missile  and  by  the  sun 
on  the  earth,  and 


AgM  =  — GM 


M 


f  1  l  ' 

~W2~  Ump  -  -y~ Ume 
V^MP  ^ME 


(5-8) 


is  the  difference  between  gravitional  forces  per  unit  mass  by  the  moon  on  the  missile  and  that  by 
the  moon  on  the  earth. 

For  a  missile  about  1,000  km  above  the  surface  of  the  earth: 

AS, 


is  less  than  one  part  per  10  millions  in  magnitude. 


5-2 


NSWCDD/MP-95/87 


A§m 

Se 


is  less  than  two  parts  per  10  millions  in  magnitude, 


where  gE  = 


is  the  magnitude  of  the  earth’s  gravitational  force  on  the  missile  per  unit  mass. 


These  values  are  at  least  one  order  of  magnitude  smaller  than  the  accuracy  of  the  best  state- 
of-the-art  accelerometer  even  in  the  laboratory  environment. 

Although  we  have  not  considered  the  effect  of  the  rest  of  the  universe  on  the  missile,  it 
turns  out  that  its  magnitude  is  much  smaller  than  the  effects  from  the  sun  and  moon. 

For  these  reasons,  (6)  is  approximated  by 


Pi2Rep  =  -^U 


R2 


Hep 


EP 


XNG 

m 


(9) 


with  negligible  error  represented  by  (7)  and  (8). 


Equation  (9)  is  the  exact  form  of  Newton’s  second  law  £F  =  ma  or  a= -  as  if  the 

m 

origin  of  the  I-frame  were  located  at  the  center  of  the  earth. 

Qualitatively  speaking,  the  magnitude  of  Ags  as  a  proportion  in  P2REp  is  neglible  because 

the  sun  pulls  both  the  earth  and  missile  together  so  that  the  relative  position  of  the  missile  with 
respect  to  the  earth  is  relatively  unaffected.  Similar  reasoning  applies  to  the  case  of  Agn,  being 
negligible. 
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6.0  ANGULAR  ROTATIONS 

In  the  addition  of  two  vectors  Vi  and  V2,  Vi  +  V2  =  V2  +  Vi  as  shown  in  Figure  6-1, 
(i.e.,  Vi  followed  by  V2  has  the  same  result  as  V2  followed  by  Vi). 


FIGURE  6-1.  VECTOR  ADDITION 

In  angular  rotations,  a  rotation  about  the  x-axes  followed  by  another  rotation  about  the 
displaced  y-axes  does  not  result  in  the  same  orientation  of  the  frame  as  a  rotation  about  the  y-axes 
followed  by  another  rotation  about  the  displaced  x-axes.  In  this  sense,  angular  rotations  do  not 
commute,  (i.e.,  do  not  act  like  vectors). 

To  see  this  by  way  of  two  simple  demonstrations,  consider  an  orthogonal  (Frame-0)  with 
Xo,  Yo,  and  Zo  axes,  as  shown  in  Figure  6-2. 

Frame  0 


Zo  Zo 


FIGURE  6-2.  O-FRAME 
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In  the  first  demonstration,  rotate  the  above  frame  about  Xo  -  axes  by  90°  (using  the  right- 
hand  rule).  This  results  in  Frame  1  (with  Xi,  Yi,  and  Zi  axes)  as  shown  in  Figure  6-3. 

Frame  1 


Yi 


Yi 

A 


Zi  *<■ 


Xi  =  X« 


FIGURE  6-3.  FRAME  1 


Next,  rotate  Frame- 1  about  Y  i  -  axes  by  90° — resulting  in  Frame-2  (with  X2,  Y 2,  Z2  axes) 
as  shown  in  Figure  6-4. 


Frame  2 


Y2=Yi  Yj  =  Yi 


FIGURE  6-4.  FRAME  2 

Now  in  the  second  demonstration,  returning  to  Figure  6-2,  rotate  Frame  0  about  Yo-axes 
(instead  of  Xo  axes)  by  90° — resulting  in  Frame  V  (with  Xi',  Yi',  and  Z\  axes)  as  shown  in 
Figure  6-5. 
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Frame  V 


FIGURE  6-5.  FRAME  V 

Next,  rotate  Frame  V  about  X  i  "-axes  by  90° — resulting  in  Frame  2'  (with  X2',  Y2',  and 
Z2'  axes)  as  shown  in  Figure  6-6. 


Frame  V 


X'2  =  X'l  X'2  =  X'l 


FIGURE  6-6.  FRAME  2' 

Obviously,  Frame  2'  (Figure  6-6)  has  a  different  orientation  from  Frame  2  (Figure  6-4).  It 
means  that,  in  describing  the  effects  of  angular  rotations,  it  is  necessary  to  specify  the  order  in 
which  a  sequence  of  rotations  is  executed. 

Consider  Frames  A,  B,  C,  and  a  vector  R,  which  has  components  RA  =  ( xA,  YA,  ZA  )T  in 
Frame  A,  RB  =  (  Xb,  Yb,  Zb  )t  in  Frame  B,  and  Rc  =  ( Xc,  Yc,  Zc  )T  in  Frame  C. 

Then 


Rb  =  CBRA 

(6-1) 

Rc  =  CgR8 

(6-2) 

Substituting  (6-1)  into  (6-2)  for  RB 
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Re=CgC°RA  (6-3) 

Since 

Rc  =  C*Ra  .  (6-4) 

From  (6-3)  and  (6-4),  we  conclude 

C£C®=C‘.  (6-5) 

Referring  to  (6-5)  and  based  on  our  discussions,  CgC®  *  C®  Cg 

As  a  matter  of  fact,  C®  Cg  has  neither  physical  nor  geometrical  meaning.  In  our  notation, 
the  second  rotation  Cg  is  written  to  the  left  of  the  first  rotation  C®  so  that  the  superscript  B  of  C® 
matches  the  subscript  B  of  Cg,  so  that  CgC®  becomes  as  in  (6-5). 

We  can  extend  this  process,  for  instance,  to: 

CEDCDcCcBC'A  =  Cl  (6-6) 

By  writing  as  in  (6-5)  or  (6-6),  we  can  check  the  consistency  of  multiplications  of  the 
coordinate  transformation  matrices.  If  superscripts  and  subscripts  do  not  match  in  a  way  we  have 
specified,  it  indicates  errors  somewhere. 
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7.0  THEOREM  OF  CORIOLIS 


Consider  two  references,  Frames  A  and  B,  with  common  origin  at  0.  Frame  B  rotates  with 
respect  to  Frame  A  with  angular  velocity  Wab-  A  vector  RqP  =  R  may  be  expressed  in  Frame  A 


and  in  Frame  B 


RA  =  IX  +  JY  +  KZ 


(7-1) 


Now,  using  (7-1) : 


R®  =  ix  +  jy  +  kz . 


,  dX  dY  dZ 
P.R  =  I— ■  + J- — l-K— 
A  dt  dt  dt 


vdl  vdJ  dK 

+  X — h  Y —  +  Z — 
dt  dt  dt 


TdX  dY  dZ 

-  I —  +  J —  +  K — 
dt  dt  dt 


(7-2) 


(7-3) 


because  the  unit  vectors  I,  J,  and  K  are  fixed  in  Frame  A,  and  do  not  vary  with  time  in  Frame  A. 
Using  (7-2): 


_  _B  . dx  .dy  .  dz  di  dj  dk 

P,R  =  i — + 1— +  k —  +  x — I-  y —  +  z —  . 
A  dt  Jdt  dt  dt  'dt  dt 


(7-4) 


Note:  In  this  case,  the  unit  vectors  i,  j,  and  k,  which  are  fixed  in  Frame  B,  rotate  relative  to 
Frame  A,  and  therefore,  are  time-varying  as  viewed  from  Frame  A. 


Now,  consider  in  (7-4).  (Referring  to  Figure  7-1). 
dt 


Since  vector  i  is  a  unit  vector,  it  cannot  change  its  magnitude.  However,  it  can  change  its 
direction  relative  to  Frame  A  because  i  is  fixed  in  Frame  B,  which  rotates  relative  to  Frame  A.  For 
an  infinitesimal  angular  displacement,  we  can  visualize  the  tip  of  the  i-vector  moving  on  the  plane 
parallel  to  the  j  -  k  plane.  Therefore,  we  may  express  Ai  caused  by  the  angular  displacement  of  i  in 
Frame  A  in  terms  of  its  displacements  in  the  j  and  k  directions. 

Referring  to  the  Figure  7-1  below,  we  see 
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W2 


FIGURE  7-1.  INCREMENTAL  ROTATION  OF  FRAME 


that  Wz  causes  Ai  =  jWzAt  in  the  j  direction  during  At  and  Wy  causes  Ai  =  -kWyAt  during  At  in  the 
-k  direction.  Summing  the  two  components,  we  have 

Ai  =  jWzAt  -  kWyAt  (7-5) 

dividing  by  At  and  taking  the  limit,  we  have 


|  =  jW,-kW, 


(7-6) 


It  turns  out  that  the  right-hand  side  of  (7-6)  is  also  equal  to  WB  x  i  as  shown  below: 


WBxi  = 


l 

W. 

1 


j  k 
Wy  W. 
0  0 


=  jWz-kWy 


(7-7) 
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It  follows: 


Similarly,  we  can  show  that: 


di 

dt 


B 


X 


i . 


JL 

dt 


=  w  Xj 


dk 

dt 


=  wB  x  k . 


(7-8) 


(7-9) 


(7-10) 


Now,  referring  to  the  right-hand  side  of  (7-4),  and  using  (7-8),  (7-9),  and  (7-10): 


di  dj  dk 

x— +  y— +  z — 
dt  dt  dt 


=  x  (wB  xi)+y  (wB  xj)  +  z(wB  xk) 

=  wB  x  (ix+jy  +  kz) 

=  wB  x  ix+ wB  x  jy-t-  wB  x  kz  (since  x,  y,  and  z  are  scalars) 

=  wB  xR® 

=  w®B  x  Rb  (since  W  =  Wab  by  definition  in  our  case)  (11) 

Now, 

PBRB  =  P  (ix+jy+kz)B 


,  dx  .  dy  ,  dz 
=  l - t-j - hk — 

dt  dt  dt  (12) 

because  i,  j,  k  are  fixed  in  Frame  B. 

So,  substituting  (7-12)  and  (7-11)  into  (7-4): 

PaRb=PbRb+W^xRb.  (13) 
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In  (13),  each  term  may  be  expressed,  after  appropriate  transformations,  in  either  Frame  A 
or  Frame  B,  (see  next  section).  So,  dropping  the  superscript  B  from  (7-13), 

PaR  =  PbR  +  WabxR  (7-14) 

which  is  called  the  Equation  of  Coriolis. 

What  the  Equation  of  Coriolis  implies  is  that  the  differentiation  (WRT  time)  of  a  vector  in 
one  frame  is  not  equal  to  the  differentiation  of  the  same  vector  in  another  frame,  which  is  rotating 
with  respect  to  the  first  frame.  And  to  obtain  the  value  of  PAR  of  (7-3)  in  terms  of  PBR  of  (7-12), 
we  have  to  add  a  correction  term  WabxR  (which  incorporates  the  effect  of  rotation  of  Frame  B 
relative  to  Frame  A)  to  PBR  as  shown  in  (7-14). 

For  some,  the  Equation  of  Coriolis  is  considered  the  second  most  fundamental  equation 
only  next  to  the  Newton’s  second  law  in  the  analysis  of  navigation  and  guidance.  Readers  who  are 
interested  in  the  geometrical  approach  in  the  derivation  of  (7-14)  may  find  it  in  other  text  books 
such  as  "Mechanics"  by  Symon  (published  by  Addison  Wesley).  Some  may  find  the  geometrical 
approach  simpler  and  easier  to  follow,  while  others  may  not.  For  this  reason,  an  analytic  approach 
is  presented  here  to  assist  the  comprehension  in  view  of  the  importance  of  the  theorem. 
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8.0  MATRIX  FORM  OF  THE  THEOREM  OF  CORIOLIS 

The  law  of  Coriolis  in  vector  form  (which  we  have  derived  previously)  between  Frames  A 
and  B  follows: 

PaR  =  PbR  +  WabxR.  (8-1) 

Coordinating  (expressing  components)  in  Frame  B: 

(PaR)b  =(PbR)b  +  W^g  xRb  .  (8-2) 

Now,  referring  to  (8-2): 

(par)‘=c;(p.r)a 

=  CJ(PRA)  because  (P,R)‘  =  (PRA)A  =  PR* 

= c;[p(c;r")] 

= cj  [  (pc*  )rb + c£  (prb  )] 

=  CB(PC*)RB  +  CBC*PRB 
=  CB(PC*)RB  +  PRB 

since  C®  C3  =C°  =  I 

Returning  to  (8-2),  recall  the  following  from  Section  1: 

W£xRb«[W“]Rb 

where  for  W®B  =  (WxWy  Wy  )T 


(8-3) 


(8-4) 
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WBK 


r  0 

-W. 

Wy 

W. 

0 

-W, 

w. 

0 

which  is  the  matrix  representation  of  the  vector  cross-product. 
Substituting  (8-3)  and  (8-4)  into  (8-2) 


CB  (PC£)RB  +  PRB  =  PRB  +  [W“]RB 


noting  that  (pbR)B  =  (PRB)B  =  PRB . 
It  follows: 


CB(PC*)RB=[W“]RB. 


(8-5) 


(8-6) 


Since  (8-6)  has  to  be  valid  for  all  RB,  we  have 

c;p(c;)=w“  (8-7) 

premultiplying  both  sides  of  (8-7)  by 

(<=:)“ =(c;)T=c; 

because  CB  ,  being  a  coordinate  transformation  matrix  (between  two  orthogonal  frames),  is  an 
orthogonal  matrix,  we  have 

C*CBP(C*)  =  C*W“  (8-8) 

since  CB  =  =  I 

PC*=C£\V“  (8-9) 

Equation  (8-9)  is  a  matrix  form  of  the  Equation  of  Coriolis. 

Since 

p/-A  _  d  pA  _  Cg(n  +  1)-Cg(n) 

B  dt  B_  At  (8-10) 
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where  CA  (n)  denotes  the  value  of  C*  at  nT. 
We  have  from  (8-9): 


CBA(n  + 1)  =  CB(n)  +  Cg(n)W^(n)At  (8-11) 

Equation  (8-10)  suggests  that  if  we  can  measure  by  gyros  fixed  to  the  body  (such  as 
the  rocket  structure),  then  we  can  update  the  orientation  of  the  body-fixed  frame  relative  to  a 
reference  frame  (such  as  the  I-frame)  by  updating  the  CA  matrix.  This  idea  is  used  in  the  so-called 

strap-down  INS  where  the  gyros  and  accelerometers  are  fixed  to  the  body  instead  of  on  the  inertial 
platform,  which  is  isolated  from  the  body  motion  by  a  servo  system  with  gimbals.  It  may  also  be 
used  in  error  analysis  of  the  EMU  platform  drift  caused  by  the  gyro-drifts  rates. 
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9.0  ROTATIONAL  FORM  OF  NEWTON’S  SECOND  LAW 


The  rotational  form  of  Newton’s  second  law  is  the  basic  starting  equation  essential  to 
describing  the  operation  of  gyro,  gyro-accelerometer,  and  base-motion  isolation  systems.  One 
such  system  is  the  stable  platform  of  the  IMU  of  a  LRBM,  which  is  controlled  by  a  servo¬ 
mechanism  using  gyros  and  gimbals. 

The  equation  also  leads  to  the  so-called  Euler  equation,  describing  the  dynamic  behavior  of 
the  rigid  body  motion,  which  is  used  in  the  SDOF  simulation  model  of  the  rocket  body  in  flight 

Consider  a  point  mass  m  located  at  point  k .  By  definition  the  moment  of  momentum  or 
“angular  momentum”  Hit  of  m  about  the  point  I  (origin  of  the  I-frame)  is 

HK=RKxmPIRK  .  (9-1) 

The  moment  of  force  or  “torque”  Mik  of  external  force  on  m  at  k  about  the  point  I  is  by  definition 
the  following: 

Mlk  =  RIkxFk  (9-2) 

Differentiating  (9-1)  WRT  time  and  noting  that  the  cross-product  of  two  identical  vectors  is 

zero: 

PiH*  =  pi(Rft  x  mPiRjk) 

=  ^R*  x  mP.R,,  +  RIt  x  mP^Rfc 
=  Rfc  x  mPjTt* 

=  R*  x  Fk 

It  follows  in  view  of  (9-2): 


PiHn,  =  (9_3) 

Summing  up  (9-3)  for  all  particles  (all  k’s),  and  defining  Hj  and  Mj  by 

H,  =  Ik  (R*  x  mkP,RIk )  (9-4) 

M,  =  Xk(Rik  xFk)  (9-5) 

and  with  the  understanding  that  Fk  refers  to  the  external  forces  only,  we  get 

P,H,  =  M,  .  (9-6) 
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This  is  the  rotational  form  of  Newton’s  second  law.  The  law  states  that  the  time  rate  of 
change  (WRT  the  1-frame)  of  the  angular  momentum  about  an  inertial  point  is  equal  to  the  torque 
about  the  same  inertial  point 

By  denoting  the  location  of  the  center  of  mass  by  C,  may  be  expressed  by  a  simple 
vector  addition: 

Rik =  Ric  +  (9-7) 

Substituting  (9-7)  into  (9-4)  and  (9-5)  and  then  the  results  into  (9-6),  after  some  algebra 
(see  for  instance  Reference  (9-1)),  we  get  the  rotational  form  of  Newton’s  second  law  about  the 
center  of  mass  (instead  of  the  inertial  point): 

P,HC  =  Mc  (9-8) 

where 

Hc=Ik(Rc>xmiP,Rcl)  (9-9) 

=  the  angular  momentum  of  the  system  of  particles  about  the  center  of  mass 
Mc  =  £kRQxFK  (9-10) 

» the  external  torque  about  the  center  of  mass. 

Equation  (9-8)  states  that  the  applied  torque  about  the  center  of  mass  is  equal  to  the  time 
rate  of  change  (WRT  the  1-frame)  of  the  angular  momentum  about  the  center  of  mass. 

The  rotational  form  of  Newton’s  second  law  is  valid  for  any  system  of  particles  whose 
internal  forces  are  central,  whether  or  not  the  body  is  rigid. 
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10.  PROTOTYPE  LINEAR  ACCELEROMETER* 


Consider  a  box  (instrument)  with  a  spring  and  a  test  mass  m.  The  box  is  rigidly  attached  to 
the  rocket  structure  as  shown  in  Figure  10-1.  The  rocket  is  moving  along  a  straight  line  under  a 
constant  acceleration. 

Under  this  steady-state  condition,  the  spring  extends  by  a  constant  displacement  x  from  P 
to  Q  along  the  unit  vector  U,  because  of  the  inertial  reaction  force  on  m,  which  is  the  reaction  to  the 
applied  force  exerted  on  P. 


u 


FIGURE  10-1.  PROTOTYPE  LINEAR  ACCELEROMETER 

The  forces  acting  on  m  at  point  Q  are  the  gravitational  forces  Fgq  (the  gravitational  force  Fq 
at  point  Q),  and  the  force  -kx  pulling  the  mass  toward  the  point  P  by  the  restoring  force  of  the 
spring. 
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Applying  Newton’s  second  law  to  m  at  Q  with  components  expressed  along  the  missile 
center  line,  by  taking  dot  product  with  u: 


mPi2RIQ,«  =  (FGQ-kx)u> 


(10-1) 


Since  by  vector  additions: 


R*  “  Rn>  +  Rfq 


(10-2) 


we  get  by  substituting  (10-2)  into  (10-1): 


(mP’R*  +  mPfRpJ •  u  =  (F^  -  kx)  •  u 


(10-3) 


Under  steady-state  conditions  Rpq  is  constant,  which  makes  PiRpq  and,  therefore,  PiRpq 
equal  to  zero.  It  follows  from  (10-3): 


PX 


u 


V  m  m  ) 


(10-4) 


Noting  that 


is  the  gravitational  acceleration  g  at  Q  (denoted  by  gQ),  we  have  from 


(10-4): 


(10-5) 


Applying  Newton’s  second  law  to  the  whole  box  considered  as  a  rigid  body  (instrument  with  the 
pivot  point  at  P ),  noting  that  the  acceleration  at  P  has  to  be  caused  by  either  gravitational  specific 
force  gp  at  P  or  nongravitational  specific  force  fNGp  at  P>  and  taking  the  component  along  the  u 
direction: 

(P,  R,p)’ u  =  (gp +fN(3>)*  u  (10-6) 

Equating  (10-5)  to  (10-6)  and  rearranging 


1"ngp  '  tJ  X  +  (gQ  gp)-u 

m 


(10-7) 
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Now  both  gp  and  gQ  are  gravitational  accelerations  (gravitational  force  per  unit  mass)  several 
thousand  miles  away  from  the  center  of  the  earth  at  points  P  and  Q  that  are  less  than  only  a  few 
centimeters  apart  Thus,  gQ  and  gp  are  essentially  equal  and  gQ  -  gp  =  0 . 

It  follows  from  (10-7) 


Equation  (10-8)  shows  that  the  mass-spring  box  (an  idealized  accelerometer)  measures  the 
externally  applied  nongravitational  force  per  unit  mass  (called  nonfield  “specific  force”),  which  is 
equal  to  the  force  (per  unit  mass)  that  the  test  mass  m  exerts  on  its  support  at  P  through  the  spring. 
That  is,  the  spring  displacement  x  is  the  measure  of  nongravitational  force  per  unit  mass  (specific 
force)  along  the  direction  u  at  point  P.  Since  x  is  proportional  to  -fNGp  •  u,  we  should  say,  more 
exactly,  that  the  accelerometer  measures  the  negative  of  the  nonfield  specific  force  along  (u). 

It  is  important  to  note  that  the  accelerometer  (a  misnomer  in  a  strict  sense)  measures 
nongravitational  specific  force  fNG,  not  f^G  +  fG ,  Newton’s  second  law  states  that 

Thus,  in  the  implementation  of  the  navigation  system  that  relies  on  the  onboard 
accelerometers,  we  have  to  always  add  the  gravitational  acceleration  (computed  generally  based  on  the 
inverse-square  gravity  model)  to  obtain  the  total  acceleration  as  shown  below  in  Figure  10-2. 


accelerometer 


pX 

total  acceleration 


fa 

G  •  computer 

FIGURE  10-2.  COMPUTATION  OF  TOTAL  ACCELERATION 

For  some,  the  fact  that  the  accelerometer  onboard  a  rocket  in  motion  under  a  gravitational 
field  measures  only  the  nongravitational  specific  force  may  seem  counter-intuitive. 

Referring  to  the  first  figure  of  this  section,  both  points  P  and  Q  located  more  than  6,000  km 
away  from  the  center  of  the  earth  are  only  less  than  a  few  centimeters  apart.  This  means  that  the 
gravitational  force  per  unit  mass  at  these  two  points  is  essentially  equal. 
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Thus,  we  can  see  that  the  effect  of  gravitational  force  on  the  displacement  of  the  spring  is 
essentially  zero.  In  another  words,  the  accelerometer  is  incapable  of  sensing  the  gravitational  force 
and  therefore  cannot  measure  it.  Now,  let  us  consider  the  combined  effect  of  thrust  and  aero¬ 
dynamic  forces  along  the  direction  u.  This  will  push  the  rocket  structure  forward  along  with  the 
pivot-point  P,  which  is  fixed  to  the  structure.  By  Newton’s  first  law,  the  mass  M  located  at  Q  will 
resist  this  change  and  try  to  remain  behind,  thus  stretching  the  spring  by  inertial  reaction.  Now, 
we  can  see  why  the  accelerometer  senses  and  measures  the  effect  of  nongravitational  force  only 
excluding  that  of  gravitational  force. 
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11.0  PULSE  INTEGRATED  PENDULOUS  ACCELEROMETER  (PIPA)* 


Consider  a  physical  pendulum  of  mass  m  with  the  moment  of  inertia  J  and  with  the  center 
of  mass  (CM)  located  at  distance  L  from  the  pivot  point  P.  It  is  damped  with  the  damping 
coefficient  C. 

The  vehicle  with  the  pendulum  (pivot  point  P)  fixed  to  it  is  under  acceleration  by  the 
combined  effect  of  both  gravitational  and  nongravitational  forces.  The  nongravitational  portion  of 
the  acceleration  is  acting  in  the  direction  labeled  as  IA  in  Figure  1 1-1. 


As  explained  in  the  case  of  linear  accelerometer,  the  gravitational  force  per  unit  mass 
or  gravitational  acceleration  at  pivot  point  is  essentially  identical  to  that  at  the  CM  of  the  pendulum. 

P 


FIGURE  11-1.  PHYSICAL  PENDULUM 

Thus,  the  gravitational  acceleration  does  not  contribute  to  the  swinging  of  the  pendulum. 
However,  the  nongravitational  specific  force  fng  (thrust  and  aerodynamic  force)  is  applied  to  the 
point  p,  which  is  fixed  to  the  vehicle,  but  not  to  CM  (unlike  the  gravitational  force  that  is  exerted  on 
both  p  and  CM).  As  p  accelerates  along  IA,  m,  according  to  Newton’s  first  law,  tries  to  remain 

stationary  in  the  inertial  space,  thus  making  angle  0  with  the  reference  line. 

Applying  the  rotational  form  of  Newton’s  second  law  about  a  pivot  axis  (into  the  paper): 


J0  =  IM 


(11-1) 


*This  section  is  adapted  from  Chapter  1.0  (the  introduction)  of  MITInstrumentation  Laboratory  Report  R-574 
(on  PIPA)  by  Kee  Soon  Chun,  March,  1967. 
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The  applied  torques  ZM  consists  of  the  following  three  parts: 

Fng 

•Torque  caused  by  the  inertial  reaction  force  -FNG=-m-^-=-mfNG  about  IA  axis. 

The  negative  sign  is  necessary  because  the  inertial  reaction  force  is  in  the  opposite 
direction  from  the  applied  Fng-  This  torque  is  from  the  diagram,  Lcos0(-mfNG)=- 
mLcos0fNG 

•Damping  (dynamic  friction)  torque  -  C0  (with  negative  sign  because  it  opposes  the 
motion). 

•Rebalancing  torque  generated  by  the  torque  generator  Mtg,  which  drives  the  angular 
deviation  0  to  zero. 

In  actual  operation,  the  angular  deviation  0  of  the  pendulum  from  the  reference  line  is 
limited  to  very  small  angles  and  constantly  rebalanced  by  the  rebalancing  torque  generated  by  the 
torque  generator  located  at  the  pivot  point  (not  shown  in  die  diagram). 

Thus,  from  (11-1),  summing  up  torques: 

J0  =  Mtg  -  mLcosOf ^  -  C0  (11-2) 

since  0  is  limited  to  very  small  angles  cos0  =  1.  Thus 

J0  +  C0  =  Mtg-mLfNG  (11-3) 

Integrating  with  respect  to  time  for  the  duration  of  At: 


At  At  At  At 


j/edt+C  J6dt=  J  M„dt-  mL  J  fNDdt 
0  0  0  0 


(11-4) 


Because  of  the  nature  of  the  rebalancing  operation,  the  angular  displacement  constantly 

reverses  its  direction.  This  causes  the  average  value  of  0  (and  average  value  of  0  as  well)  to 
approach  zero  for  the  time  duration  much  greater  than  the  period  Ttg  of  each  rebalancing  torque 
pulse,  which  is  indeed  the  case  in  actual  application. 

Thus,  for  At  >  >  T, 


J*At .  f  At .. 

0dt=  I  0dt  =  O 
0  Jo 


(11-5) 
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It  follows  from  (1 1-4)  with  (1 1-5): 


Jo  —  J0  ^NG^t  • 


(11-6) 


In  the  actual  system,  Mtg  is  proportional  to  the  constant  DC  current  I,  through  the  torque 
generator  (with  the  permanent  magnet).  Thus, 


Mtg  —  Stg  I 


(11-7) 


where  Stg  is  the  torque  generator  sensitivity.  It  follows, 


(11-8) 


By  making  the  integration  duration  At  sufficiently  small  ( while  keeping  At »  Ttg )  relative 
to  the  thrust  and  aero-dynamic  force  variation,  we  may  treat  fng  constant  during  At  Thus, 


dt  =  fNG 


At  =  AVng 


(11-9) 


where  AVng  is  the  velocity  increment  caused  by  fng  during  At .  Substituting  (11-8)  and  (11-9)  into 
(11-6): 


StgIAt  —  mLAVNG 


or 


AVK0=-^IAt. 

mL 


(11-0) 


Thus,  AVng  is  proportional  to  I  At  (integration  of  current),  or  the  amount  of  electric  charge 
passed  through  the  DC  torque  motor  during  At .  If  we  use  electric  pulses 

with  constant  magnitude  (instead  of  DC  current),  the  summation  (integration)  of  pulses  is  also 
equal  to  the  amount  of  electric  charges.  Because  AVNG  is  the  integration  of  fNG  (which,  in  turn, 
proportional  to  the  summation  (integration)  of  current  pulses)  during  At,  this  type  of  instrument  is 
called  pulse  integrated  pendulous  accelerometer  (PIPA). 
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From  (11-10)  with  At  fixed 


or 


d(AVj  =  ^<U 

_  A VNG  „  (using  (10)  again)  (11-11) 


d(AVN0)  dl 

AVno  I  (11-12) 


Equation  (11-12)  shows  that  one  part  per  million  error  in  the  DC  current  measurement 
would  cause  the  corresponding  one  part  per  million  error  in  the  estimation  of  the  velocity  increment 
AVng.  In  actual  implementation,  the  constant  current  source  for  the  DC  motor  is  achieved,  in 
essence,  by  applying  a  constant  voltage  from  a  zener  diode  across  a  precision  resister.  But,  both 
zener  diodes  and  precision  resistors  are  not  only  temperature  sensitive  but  also  deteriorate  with 
time.  For  this  reason,  to  achieve  more  accuracy,  PIPA  is  replaced  by  PIGA  in  the  newer  long- 
range  ballistic  missile.  PIGA  is  discussed  in  a  later  section. 
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12.0  ROTATIONAL  MOTION  OF  THE  MISSILE  BODY 


In  the  previous  section  on  the  rotational  term  of  Newton’s  second  law  (Section  9.0),  the 
following  equation  has  been  derived: 

PiHc  =  Mc  (12-1) 

where  Hc  is  the  angular  momentum  about  the  center  of  mass,  Me  is  the  applied  moment  (torque) 
about  the  center  of  mass,  and  Pi(-)  is  “£■(')  in  the  1-frame. 

Applying  the  theorem  of  Coriolis  to  (12-1)  between  the  I-frame  and  a  body  fixed  body  (B)- 
ffame  of  the  missile  and  expressing  the  components  in  the  B-frame: 

(P,HC)B  =  P,H?  +  W|  x  H?  =  M?  (12-2) 

As  shown  in  Section  9.0,  the  angular  momentum  He  is  given  by 

Hc  =  XmkRCK  x  PjRck  (12-3) 

Applying  the  theorem  of  Coriolis  to  Rck  between  the  I-frame  and  the  B-frame,  which 
rotates  relative  to  the  I-frame  with  the  angular  velocity  Wib,  we  have: 

P[Rck  =  PbRck  +  Wib  x  Rck 

-  Wib  x  Rck 

since  PbRck  =  0  because  Rck  is  fixed  in  the  B-frame.  (Strictly  speaking,  rocket  exhaust  gases 
cause  the  center  of  mass  of  the  rocket  to  shift,  and  Rck  is  not  fixed  in  the  B-frame.  The 
derivation  here  assumes  that  the  shift  of  the  center  of  mass  is  negligible  for  the  time  span  of  this 
discussion.) 

Substituting  Wib  x  Rck  for  PiRck  in  (12-3),  we  have 

(12-3b) 

K 


Using  a  vector  identity: 

V,x(V2xV3)=(V1V3)Vj-(V1V2)Vj 

=  (V,  'V1)V2-(V1V!)V1  for  V,=V3 
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we  get  from(12-3b): 

Hc  =  Xmk(^CK  '  RckJWjb  —  ^mk(R-CK  ’  Wm)RcK  • 


(12-4) 


At  this  point,  we  define  the  x,  y,  z  axis  of  the  B-frame  with  the  origin  at  center  of  gravity  (CG)  of 
the  missile  as  shown  in  Figure  12-1. 


FIGURE  12-1.  B-FRAME  FIXED  TO  THE  MISSILE 
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where  the  axes  are  defined  as : 

X  axis —  Along  the  missile  center  line 
Y  axis  — Perpendicular  to  the  center  line  at  the  CG 

Z  axis  — Perpendicular  to  both  the  x  and  y  axis  at  the  origin  so  as  to  form  an  orthogonal 
triad  by  the  right-hand  rule  (X  ®  Y  =  Z) 

Now  consider  the  transversal  cross-section  of  the  missile  on  the  y-z  plane  as  shown  below  in 
Figure  12-2. 


z 


FIGURE  12-2.  CROSS-SECTION  (y-z  PLANE)  OF  THE  MISSILE 
We  see  that  for  any  point  Pi(yi,zi),  there  exist  P2(y2,Z2)>within  the  missile  such  that 

yjZ,  +  y2z2  =  ab  +  (-ab)  =  0 . 

Extending  this  situation  to  cover  the  entire  volume  of  the  missile, 
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Next,  consider  a  longitudinal  cross-section  on  the  x-y  plane  as  shown  below  in  Figure  12-3. 


x 


We  see  that  for  any  point  P3(x3, 73),  there  exist  P^x^  y4),  within  the  missile  such  that 

x3y3 + x4y4 = cd + (~cd) = 0. 


Therefore,  for  the  entire  volume  of  the  missile, 

Xx*yK=o- 

K 


(12-6) 


Similarly  on  the  xz  plane: 


5^xkzk  —  ® • 

K 


(12-7) 


By  denoting  and  W®  in  (12-4)  by 
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Y  > 
XK 

fw^ 

50 

fi¬ 

ll 

yK 

and  = 

Wy 

^ZK> 

IwJ 

we  have 


!>,  -RkK-Im,  (x' +yi+4) 

k 


'W:  ' 

w 


,w  , 

V  z ) 


Using  (12-5),  (12-6),  and  (12-7)  in  the  second  term  of  (12-4): 


Xmk(Rqc  •W.bJRck  —  mk  (xKWx  +  yKWY  +  zKWz) 


yK 


rIm.x|W; 

2>*y|W, 


Substituting  (12-9)  and  (12-10)  into  (12-4): 


fWl 

H 1  = 

£mk(x^  +  z^)Wy 

= 

JyWy 

£mk(4+y^)Wz; 

Iw 

which  may  be  written  as  : 


O 

o 

K 

fw‘l 

H?  = 

o 

<— 1 

o 

Wy 

N 
H- > 

o 

o 

IwJ 

«— 1 

H 

II 

M 

1 

(y|+4) 

(h 

0 

0^ 

Jy=£mK| 

(x£  +  Zk)  and  J  = 

0 

Jy 

0 

Jz  =  ^mK( 

xK  +  yx) 

lo 

0 

Jz; 

(12-13) 
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In  (12-12)  and  (12-13),  we  note  that  all  off-diagonal  elements  of  the  J-matrix,  which  is 
called  the  moment  of  inertia  tensor,  are  all  zeros,  that  is,  the  J-matrix  is  a  diagonal  matrix,  thus 
greatly  simplifying  analysis  and  computations.  The  choice  of  body  axes,  which  result  in  the 
diagonal  J-matrix  is  called  the  "principal  axes."  What  seems  to  be  a  natural  choice  of  axes  in  our 
case  happens  to  coincide  with  the  principle  axes. 

For  analytical  methods  (which  are  not  simple)  for  determining  the  principal  axes,  refer  to 
Mechanics  by  K.R.  Symon  or  Classical  Mechanics  by  H.  Goldstein. 

Two  proven  theorems  for  selecting  the  principal  axes  follow: 

•  Any  axis  of  symmetry  of  a  body  is  a  principal  axis  such  as  the  x-axis  in  our  case. 

•  Any  plane  of  symmetry  of  a  body  is  perpendicular  to  a  principal  axis.  In  our  case,  the 
plane  (containing  the  x-axis),  which  is  perpendicular  to  the  y-axis  is  perpendicular  to  the 
principal  axis,  thus  making  the  y-axis  the  principal  axis.  Similarly  for  the  z-axis. 

Now  referring  to  (12-2),  we  have,  using  (12-11): 


PbH 


JyWy 

= 

JyWy 

UwJ 

Iwj 

(12-14) 


since  Jx,  Jy,  Jz  are  constants  in  view  of  (12-13)  for  constant  fixed  mass  missiles.  (Jx,  Jy,  Jz  are 
not  constants  for  variable  mass  missiles  such  as  real-world  missiles  in  flight) 

And, 


W“xH“  = 


l 

Wx 

IW 


w 

»Ty 

J  w 

Jy  y 


k 

W, 

W 


'JzWywz-jywywz" 

'(J.-J,)w,w.' 

JZWZWZ-JZWZWZ 

= 

(J„-J,)w,w, 

(12-15) 


Substituting  (12-14)  and  (12-15)  into  (12-2),  and  expressing  the  components  of  M®  by  Mx,  My, 
Mz  we  have 
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(12- 16a) 

(12-16b) 

(12-16c) 

Equation  (12-16)  is  called  Euler’s  equations  (of  motion)  for  a  rigid  body  (along  the 
principal  axes). 

Equation  (12-16)  may  be  solved  for  Wx(t),  Wy(t),  Wz(t),  (which  are  the  angular  velocities 
of  the  missile  body  relative  to  the  inertial  space)  in  terms  of  Mx,  My,  Mz.  Mx,  My,  and  Mz  are  the 
externally  applied  torques  generated  by  the  reaction  force  to  thrust  (which  is,  by  custom,  simply 
called  thrust)  and  the  aerodynamic  force  about  the  CG  of  the  missile  (known  as  a  function  of  time 
in  simulation). 

In  Section  17.0  on  the  matrix  form  of  the  theorem  of  Coriolis,  we  have  derived 


jLr*1  —  p1  wBK 

^  V'B  ~  '-'B  ”IB 


(12-17) 


(with  the  superscript  A  of  (12-9)  in  that  section  replaced  by  I)  where  W®K  is  a  3X3  matrix  given  by 


W“  = 


O 

w, 

-w„ 


-w. 

Wy 

0 

-W, 

Wx 

0 

(12-18) 


With  Wx(t),  Wy(t),  Wz(t)  from  (12-18),  then  Equation  (12-17)  may  be  solved  for  ClB  (t) 


cL(0= 


Qi(0  Q2W  QaM 

^21  (t)  ^22(1)  C23  W 
C3,(t)  C32(t)  C33(t) 


(12-19) 


which  is  the  coordinate  transformation  matrix  from  the  B-frame  to  the  I-frame,  thus  specifying  the 
angular  orientation  of  the  missile  in  the  inertial  space. 


Newton’s  second  law  applied  to  the  CG  of  the  missile  and  expressed  component-wise  in 
the  inertial  frame  is: 


Pi  Xj  —  fGX  +  fNGX 

P’Y^fev+fN 


ANGY 


(12- 20a) 
(12-20b) 
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—  foz  +  ^NGZ 


(12-20c) 


where  P*  ( •) = — -  ( •)  in  the  inertial  space,  fg  is  the  gravitational  force  per  unit  mass,  and  f^G  is 


dt 


the  nongravitational  specific  (per  unit  mass)  force  which  includes  both  thrust  and  aerodynamic 
forces.  The  fNG  is  most  conveniently  expressed  in  the  B-frame  and  transformed  to  the  inertial 


frame. 


The  three  equations  of  (12-20)  completely  specify  the  translational  motion  (position  and 
velocity)  of  the  missile’s  CG  along  the  inertial  Xi,  Yi,  Zi  axes,  while  another  three  equations 
(12—16)  completely  specify  the  angular  orientation  (which  is  the  integration  of  the  angular  rates)  of 
the  missile  in  the  inertial  space.  For  this  reason,  the  mathematical  model  for  simulation  based  on 
these  six  equations  are  called  a  six  degree-of-freedom  (6  DOF)  model. 

If  we  are  concerned  only  with  the  translational  motion  of  the  CG,  ignoring  the  missile 
orientation,  only  the  three  equations  of  20 ’s  are  used,  and  its  model  is  called  a  3  DOF  model. 

Since  the  motion  of  CG  of  the  missile  is  the  same  as  the  motion  of  the  point  mass  in  which 
the  entire  mass  of  the  missile  is  concentrated  at  the  CG,  it  is  also  called  the  point  mass  model. 
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13.0  THE  GYROS* 


The  rotational  form  of  Newton’s  second  law  about  the  center  of  mass  C  is,  as  derived 
previously : 


P,HC  =  MC. 


(13-1) 


FIGURE  13-1.  GYRO  ELEMENT 

The  assembly  (the  spinning  rotor  and  its  support  gimbal)  shown  in  Figure  13-1  is  called 
“gyro  element”  (ge),  with  axis  of  the  ge  frame  as  shown  in  the  figure.  Replacing  subscript  C  by  ge 
(cm  of  ge)  in  (13-1)  and,  applying  the  law  of  Coriolis  between  I-frame  and  ge  frame, 

PlHge  =  PgeHge  +  Wj.ge  x  Hge  =  Mge  (13-2) 

where  the  subscript  i-ge  is  the  same  as  I-ge. 

The  angular  momentum  of  the  gyro  element  Hge  consists  of  the  spin  part  Hs  (constant)  and 
nonspin  part  H„s,  that  is  : 


Hge=Hs  +  Hns  (13-3) 

where  Hm  is  the  angular  momentum  vector  of  the  gyro  element  other  than  that  of  the  spin  rotor. 
Substituting  (13-3)  into  (13-2): 
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PgeHs  +  PgeHns  +  Wi.ge  X  (Hs  +  Hns)  =  Mge  .  (13-4) 

In  (13-4),  PgeHs  =  0  because  Hj  is  constant  and  fixed  in  ge  frame,  and  therefore 

Wi.gex(Hs  +  Hns)  =  WiexHs 


because 

Hg  »  Hns  • 

It  follows  that 

PgeHns  +  Wj.ge  X  Hs  =  Mge  • 

(13-5) 

When  PgeHns  is  negligible  compared  to  Wi.ge  x  Hs,  (13-5)  reduces  to 

Wi.ge  X  Hs  ”  Mge 

(13-6) 

Equation  (13-6)  is  the  approximate  performance  equation  for  the  practical  gyro.  The 
geometric  interpretation  of  (13-6)  is  depicted  in  Figure  13-2  with  axes  defined  in  Figure  13-3. 
That  is,  the  spin  angular  momentum  Hg  precesses  (rotates)  relative  to  the  inertial  space  with  the 
angular  velocity  of  ge  (relative  to  I-frame),  Wi.ge,  in  an  attempt  to  align  itself  with  the  applied 
torque  vector  Mge  about  the  center  of  mass  of  ge. 
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GYRO 

PRECESSION 

VECTOR 

wi(ge) 


^i(ge) x  Hs  -  M 


The  spin  momentum  of  a  gyro  precesses 
relative  to  inertial  space  in  an  attempt 
to  align  itself  with  die  applied  torque. 


''applied 


APPLIED  FORCE 
PRODUCING  TORQUE 


TORQUE 

VECTOR 

M 


GYRO 

PRECESSION 


GYRO  SPIN  ANGULAR 
MOMENTUM  VECTOR 


Hs 


FIGURE  13-2.  REPRESENTATION  OF  THE  BASIC  LAW  OF  MOTION 


FIGURE  13-3.  DEFINITION  OF  THE  PRACTICAL  GYROSCOPE  AXES 
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Referring  to  (13-5),  we  express  each  vector  by  its  components  by 


'WxV*e 


Wge  = 

i-ge 


w 


w  , 

\  Z  / 


H  = 


'0  ' 

0 

vH.y 


and 


P  H  =P 

ge  ns  ge 


(J  W  ) 

gc 

(J  P  W  ) 

gc 

Jx  P  Wx' 

X  X 

J  w 

x  gc  X 

J  P  wv 

J  P  w 

y  y 

y  y 

k  wj 

y  g«  y 

k  pg.  wj 

J  P  Wi 

\  z  z  J 

M  = 

ge 


(Mx)*e 
M 

y 

.M  . 

V  Z  ) 


where  Jx,  Jy,  and  Jz  are  the  moments  of  inertia  of  gyro-element  about  its  x,  y,  and  z  axes.  (J  is 
used  instead  of  I,  which  is  used  for  identity  matrix.)  Substituting  these  into  (13-5),  and 
performing  the  cross-product  Wi .  ge  x  Hs,  we  get  the  following: 


JxPWx  +  HsWy  =  Mx 
JypWy-HsWx  =  My 
JzpW2  +  0  =  Mz . 


(13-7a) 

(13-7b) 

(13-7c) 


Understanding  the  components  are  expressed  in  ge-frame. 

Of  these,  only  13-7b  is  of  interest  to  us,  which  is  repeated  in  a  modified  form  with  original 
subscripts  for  Wx  and  Wy. 


P(JyW(i.ge)y)=  My  +  HSW (i-ge)x  •  (13-8) 

Equation  13-8  shows  that  the  time  rate  of  change  of  the  angular  momentum  about  the  y-axis 
of  gyro-element  is  equal  to  the  sum  of  the  applied  torque  about  the  y-axis  and  the  gyro-scopic 
torque  created  by  the  interaction  of  the  gyro’s  spin  angular  momentum  Hs  and  the  angular  velocity 
of  the  gyro-element  about  its  x-axis. 

Some  authors  refer  to  this  gyroscope  torque  as  the  “reaction”  torque  caused  by  the 
Wfl.ge)X.  However,  it  is  not  “reaction”  in  the  sense  of  Newton’s  third  law  because  HsW(i.ge)X 
acts  in  the  y-direction,  while  W(i.ge)X  is  in  the  x-direction,  which  is  perpendicular  to  the  y- 
direction. 
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FIGURE  14-2.  THE  GYROELEMENT 


XTg“  INPUT  AXIS 
IA 


FIGURE  14-3.  GYRO  ELEMENT  COORDINATE  AXES 

Xgu  along  IA  (input  axis)  of  gyro-unit  (case) 

Ygu  along  OA  (output/precession)  axis  of  gu  (case) 
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Zgu  along  spin  reference  axis  (SRA)  of  gu  (case) 

When  a  =  0,  Zgu  =  Zge  or  SRA  =  SA. 

a  is  the  angle  the  spin  axis  (SA)  makes  WRS  the  SRA,  which  coincides  with  the  Zgu)  as  a  result  of 
the  rotation  of  the  gyro-element  about  OA. 

Now 

Wj.ge  =  WLgu  +  Wgu.ge  .  (14-1) 

Expressing  (14-1)  in  the  ge  frame: 

W5.-CJW' +w;, .  a4.2) 

Referring  to  Figure  14-3  and  noting  that  for  a  small  angle  a,  sina  =  a  and  cosa  =  1, 


Xgg  —  Xgu  cosoc  +  Zgu  since  —  Xgu  +  ocZgu 


Yge  =  Ygu 


Zge  — Xgusina  +  Zgucosa  —  — ’Xg^cx  Zgy 


(14- 3a) 
(14-3b) 

(14- 3c) 


From  the  above  equations, 


( Y  V  / 


v  ZK  . 


1  0  a 
0  1  0 
-a  0  1 


Ve/V  Y° 


g» 


go 


go 


V^go  J 


(14-4) 


Thus,  C®*  in  (14-2)  is  given  by 


C£  = 


f  1  0 
0  1  0 


-a  0  1 


(14-5) 


We  designate  the  components  of  Wfege  and  W® “u  as  follows 


W JL  = 


Wy 


(14-6) 
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W&  = 


wOA 

vWSRAy 


(14-7) 


Now,  ge  can  rotate  with  respect  to  gu  about  Yge  =  Ygu  =  OA  only,  creating  the  angle  a. 

'  0  ' 


wge  = 

gage 


Pa 

X®  , 


(14-8) 


Substituting  (14-5),  (14-6),  (14-7),  and  (14-8)  into  (14-2) : 


Wx  =  WiA  +  aWSRA 

(14-9) 

Wy  =  W0A  +  Pa 

(14-10) 

Wz  =  -aWiA  +  WSRA  • 

(14-11) 

The  applied  torque  My  about  the  OA  ( Ygu  =  Yge  axis ),  consists  of : 


-Ka — Elastic  torque  proportional  to  angular  displacement  a  and  opposing  a 


-CPa — Damping  (  dynamic-friction  )  torque  proportional  to  angular  velocity  Pa,  and 
opposing  Pa 

Mtg — Commanded  torque  applied  by  the  torque  generator  located  on  the  output  axis 


Thus,  we  may  write 


My  =  -Ka  -  CPa  + 


(14-12) 


By  designer’s  choice,  when  Ka  predominates  (  relative  to  CPa  term  ),  the  instrument 

becomes  a  rate  gyro.  If  the  CPa  term  predominates,  the  instrument  becomes  a  (rate)  integrating 
gyro.  This  is  explained  below: 

In  Section  12.0,  we  derived  the  following  equation: 

JyPWy  -  HsWX  =  My  (14-13) 

substituting  (14-9),  (14-10),  and  (14-12)  into  (14-13): 


JyP(WOA  +  Pa)  -  Hs(W!a  +  aW SRA)  =  -Ka  -  CPa  +  Mtg  (14-14) 
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which  gives  (expressed  in  the  ge  frame): 


JyP2a  -  HsWu  +Mtg  -  Ka  -  CPa  +  HsaWSRA  -  JyPW0A  •  (14-15) 

For  the  operational  environment  in  which  SDOF  gyros  are  used,  (ocWsra  )  and  PWoa  are 
negligible  and,  therefore,  the  last  two  terms  on  the  RHS  of  (14-15)  may  be  ignored.  It  follows 
from  (14-15 ) : 

JyP2a  =  HsWia  -  Ka  -  CPa  +  Mtg.  (14-16) 

Equation  (14-16)  shows  that,  in  the  absence  of  the  commanded  torque  Mtg,  the  gyroscopic 
torque  HsWia  causes  angular  acceleration  P2a,  which  overcomes  the  elastic  torque  Ka  and  the 
damping  torque  CPa. 

14.1  RATE  GYROS 

For  rate  gyro  applications,  Mtg  =  0.  From  (14-16): 

Jy  ,  C  H. 

l^P2a+lc  Pa+a=~Kf’WiA  (14-17) 

By  design  K»C»Jy  for  a  rate  gyro.  So,  for  a  constant  input  value  of  Wja,  a  settles 
rather  quickly  to  the  final  value  of 


a  = 


AW1 

K 


IA 


(14-18) 


J  C 

because  <  <  —  <  <  1  in  (17). 


Since  the  output  angle  a  is  proportional  to  the  angular  rate  Wia  of  the  gyro  unit  (case) 
relative  to  the  inertial  space,  the  SDOF  gyro  of  this  type  is  called  “rate  gyro.” 


14.2  RATE  INTEGRATING  GYRO  (RIG) 

For  a  RIG,  the  elastic  constant  is  removed,  thus  making  K  =  0,  and  the  damping  constant 
C  is  made  large  (larger  than  that  of  the  rate  gyro). 

Thus,  from  (14-16)  setting  K  =  0, 

JyP2a  +  CPa  =HsWia  +  Mtg  (14-19) 

or 
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=  — WIA  +-Mt, 

^  1A  y-,  ti 


(14-20) 


For  a  constant  input  value  of  WiA,  the  —  Pa  term  vanishes  rather  quickly  (relative  to  a) 
J 

because  G»Jy  by  design  or  -pf-  <  <  1 ,  and  (14-20)  reduces  to  (in  the  absence  of  Mtg), 


Pot  =  --i-W 

m  c  u '  (14-21) 

Since  P  ( •) =-^  ( •),  we  have  from  (14-21) 


(14-22) 


According  to  (14-22),  the  angular  displacement  a  of  ge  relative  to  gu  (case)  about  the 
output  axis  is  proportional  to  the  integration  with  respect  to  time  of  the  angular  rate  Wu.  Hence, 
the  SDOF  gyro  of  this  type  is  called  a  “rate  integrating  gyro”  or  simply  an  “integrating  gyro.” 


Denoting  J  Wu  dx  =  Ww  t=  0  u  in  (14-22)  for  a  constant  Wia,  we  have 
o 


a  =  a 


IA 


It  is  to  be  noted  that  ccqa  is  not  exactly  equal  to  Oja,  but  scaled  by 


(14-23) 


14.2.1  Gvro  Transfer  Function 

(Those  who  are  not  familiar  with  the  Laplace  Transform  method  may  skip  this  section) 

By  Laplace  Transform,  differential  equations  are  transformed  to  algebraic  equations.  This 
sometimes  facilitates  the  system  analysis,  without  computer  simulations  of  the  differential  equation 
in  the  time  domain. 

Denoting  the  Laplace  Transform  L  of  a(t)  by  A(s),  that  is,  L[a(t)]  =  A(s)  and  recalling 
that  L[Pa(t)]  =  SA(s)  assuming  a  (0+  )= 0 

and 


14-6 


NSWCDD/MP-95/87 


L[P2a(t)]  =  s2A(s)  assuming  a(0+)=0  Pa(0+)=0,  we  get  from  (14-16),  with  Mtg  =  0, 

Jys2a(s)  +  Csa(s)  +  Ka(s)  =  HsWia(s)  =  (JyS2  +  Cs  +  K)a(s)  =  HsWia(s)  (14-24) 
which  gives 


a(s)= 


H 


Jys  +Cs+K 


WIA(s) 


(14-25) 


or 

a(s)  =  G(s)Wia(s)  (14-26) 

where 


G 


(s)= 


<*oa  (»)_  H 
Wia(s)  Jys2+Cs+K 


(14-27) 


G(s)  is  called  the  transfer  function,  which  relates  the  input  Wia(s)  to  the  output  a(s)  =  ccoa(s) 
with  subscript  OA  indicating  that  a  is  the  rotation  about  the  output  axis. 

For  a  constant  step  input  of  Wia,  Wu  ( s ) = — . 

s 

Applying  the  final  value  theorem  of  lim  a  ( t ) = lim  sa  ( s ) ,  we  have,  from  (14-25), 

t-oo  *  -0 


a  (t)  =  lims 

t-oo  s-0 


H, 

Jys2+Cs+K 


W 


IA 


(14-28) 


or 


a(t)  =  ^-Wu  as  t  ->  oo. 

K  (14-29) 

which  is  identical  with  (14-18)  for  the  rate  gyro.  Recalling  that  K»C»Jy  for  a  rate  gyro  (that  is, 

C  is  negligible  relative  to  K,  yet  C/Jy  is  big),  we  can  see  from  (14-28)  that  a(t)  approaches  the 
final  value  faster  as  S  approaches  zero  than  the  case  when  C  is  comparable  to  K  or  bigger  than  K. 

For  the  integrating  gyro  K  =  9,  while  keeping  C»Jy.  For  this  case,  from  (14-25) 


a(s)  = 


(14-30) 
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The  transfer  function  in  (14-30)  has  a  pole  at  the  origin  (s  =  0),  and  therefore,  we  cannot 
apply  the  final  value  theorem. 

For  a  constant  step  input  of  Wu ,  Ww  (s)=y  .  So,  from  ( 14-30) 


where 


F(s)  = 


H.  W1A 
(JS+C)  s 


(14-32) 


By  denoting  the  inverse  Laplace  Transform  of  F(s)  by  L-iF(s)  =  f(t) 


(14-31) 


f(t)  =  L-1 


(14-33) 


It  follows,  using  the  mathematical  table  for  Laplace  Transform : 


L  J 


(14-34) 


In  a  RIG  C»J  by  design  or  y-  >  >  1 ,  so  we  may  regard  exp  (-  j  t )  to  approach  zero  very 
quickly.  It  follows  from  (14-34)  for  a  step  input  of  Wl\: 


(14-35) 


Using  a  Laplace  Transform  identity  in  (14-31)  with  (14-35), 
o(t)  =  L-‘o(S)  =  L-‘|F(s) 
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,  H  WIA  H  ,  H 

=  J  ~^dT=-^ Wu  j dT~ Wut 


a(t)=^-eIA 


(14-36) 


where  0^  is  the  angular  rotation  about  IA  caused  by  Ww  for  the  duration  of  t,  or  0^  =  WM  t. 
Note  that  (14-36)  is  identical  with  (14-23),  for  the  RIG. 
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15.0  STABLE  PLATFORM 
(IMU  PLATFORM) 


The  stable  platform  of  an  LRBM  is  a  platform,  which  is  kept  aligned  with  the  earth- 
centered  1-frames  by  a  gimballed  servo-mechanical  control  system,  which  uses  gyros  to  detect 
rotational  motion  of  the  stable  platform.  Three  mutually  orthogonal  accelerometers  are  rigidly 
attached  to  the  stable  platform  and  are  used  to  determine  the  position  and  velocity  of  the  missile. 
The  stable  platform  with  the  accelerometers  and  the  associated  electronics  is  called  an  inertial 
measurement  unit  IMU  and  is  the  central  part  of  the  missile  guidance  system. 

First,  we  explain  the  single-axis  stable  platform  before  discussing  the  more  general  three- 
axes  stable  platform.  Consider  a  conceptual  single-axis  stable  platform  (plate)  with  the  body-fixed 
axes  as  shown  in  Figure  15-1. 


FIGURE  15-1.  CONCEPTUAL  SINGLE-AXIS  STABLE  PLATFORM 

The  Xg-axis  is  constrained  to  be  parallel  to  the  Xi  axis  of  the  I-ffame.  The  plate,  which 
is  perpendicular  to  the  Xg-axis,  is  subject  to  the  rotational  disturbances  about  the  Xg-axis  caused 
by  unwanted  disturbances.  This  means  that,  although  the  plate  itself  has  a  fixed  orientation  in  the 
inertial  space,  the  Yg-axis  and  Zg-axis  fixed  on  the  plate  may  deviate  from  the  original  Yi  -  axis 
andZi  -  axis. 
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The  control  problem  is  how  to  counterbalance  the  deviation  of  the  Yg-axis  and  Zg-axis 
caused  by  the  unwanted  disturbance  rotation  about  the  Xg-axis,  so  that  the  Yg-axis  and  Zg-axis 
return  to  the  directions  parallel  to  the  Yj  -  axis  and  Zi  -  axis. 

This  may  be  achieved  by  a  scheme  as  shown  in  Figure  15-2.  The  center  piece  of  the 
scheme  is  an  SDOF  gyro  with  the  gyro  unit’s  IA  designated  as  the  X-axis,  the  OA  as  the  Y-axis, 
and  the  SRA  as  the  Z-axis.  SRA  is  the  direction  of  the  spin  axis  (SA)  when  SA  is  orthogonal  to 
IA.  SRA  is  the  reference  direction  of  SA.  SA  is  always  orthogonal  to  OA. 


GYRO 

ELEMENT 


SIGNAL  GENERATOR 


0, 

OUTPUT 

AXIS 

GYRO  UNIT 


FIGURE  15-2.  MECHANICAL  SCHEME  FOR  SINGLE-AXIS  STABLE  PLATFORM 

The  type  of  gyro  used  for  this  purpose  is  called  a  RIG  or  simply  an  integrating  gyro  for  the 
following  reasons: 

The  IA  of  the  gyro  unit  is  constrained  by  the  motor-shaft  structure  not  to  be  rotatable  with 
respect  to  the  base  except  by  the  servo-motor  control.  When  the  base  with  the  gyro  on  it  (with  the 

servo  control  off)  rotates  about  the  +X  axis  (IA)  by  an  angle  0X,  the  SA  will  precess  (turn) 
downward  toward  the  +X  direction,  following  the  gyroscopic  principle  that  the  spin  vector  will  try 
to  align  itself  with  the  torque  vector  (the  direction  of  rotation). 

This  precession  of  SA  rotates  the  gyro-element’s  SA  axis  by  an  angle  0y  (from  SRA 
direction)  about  the  OA  (+Y  direction). 
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For  the  RIG,  by  special  design  construction,  0y  is  related  to  0X  with  close  approximation 
for  small  angles  by:  (see  a  previous  section  on  SDFG) 


(15-1) 


in  which  Hg  is  the  constant  spin  angular  momentum  of  the  gyro  spin  rotor,  and  C  is  the  dynamic 
friction  (so  called  damping)  coefficient  of  the  gyro  element  about  the  output  axis  (OA).  Equation 

15-1  shows  that  the  output  angle  0y  about  OA  is  proportional  to  0X  or  to  the  integration  of  the 


angular  rate  0x  about  I  A  which  is  the  input  angular  rate.  For  this  reason,  the  device  is  called  a 


RIG. 

The  0y  about  OA  caused  by  0X  about  IA  is  sensed  by  the  signal  generator  (SG) 
mounted  on  the  output-axis.  This  electrical  signal  from  SG  drives  the  servo-motor  whose  shaft 
coincides  with  the  IA  of  the  gyro,  so  that  the  gyro  unit  is  rotated  in  the  opposite  direction  from  the 

initial  0X  rotation. 


This  operation  turns  the  spin  axis  (SA  of  gyro)  upward  until  0y  becomes  zero  or  0y  is 
rebalanced,  at  which  time  the  signal  from  SG  becomes  zero,  causing  the  motor  operation  about 
IA  to  cease.  With  this  control  scheme,  Yb  and  Zb  axis  misoriented  by  the  initial  rotation  0X 

about  IA  of  the  platform  is  returned  to  the  original  Yi  and  Z\  axes  by  the  counter-rotation  (-0X) 
about  the  (+/A)  direction  of  the  platform. 

Referring  to  Figure  14.2,  suppose  the  SA  drifts  downwards  generating  0y  in  the  absence 
of  rotation  about  the  input  axis  (Xs-axis),  even  though  Yb  and  Zb  axis  point  to  the  correct 
directions,  and  thus  do  not  need  correction. 

Nonetheless,  in  the  actual  system,  the  SG  will  detect  0y  and  make  the  servo-motor  turn  the 

platform  about  IA  in  such  a  way  as  to  null  the  0y.  This  would  cause  the  Yb  and  Zb  axes  to 
deviate  from  the  correct  directions  when  correction  is  not  necessary.  This  results  in  the  platform 
drift-error  caused  by  undesired  drift  of  the  gyro  spin  axis.  This  phenomenon  is  called  gyro  drift. 
By  custom,  the  angular  rate  of  gyro  drift  is  simply  called  gyrodrift 

Consider  a  hypothetical  platform  (block)  with  three  single-axis  platform  control  systems 
with  each  SDOF  gyro's  IA  ,  OA,  and  SRA  oriented  as  shown  in  Figure  15.3. 
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FIGURE  15-.3.  THREE-AXIS  STABLE  PLATFORM 
This  is  a  hypothetical  arrangement  for  tutorial  purposes. 

At  t  =  0,  we  assume  the  XB,  YB,  ZB  axes  of  the  platform  are  parallel  to  the  Xl5  Yj,  Zt  axes  of 
an  inertially  nonrotating  frame. 

The  input  axis  of  the  X-gyro  designated  by  IAX  points  to  the  inertially  fixed  direction  if 
there  are  no  rotations  about  IAy  (IA  of  Y-gyro)  and  IAZ  (IA  of  Z-gyro).  Similarly,  IAy  points  to 
the  inertially  fixed  directions  if  there  are  no  rotations  about  IAX  and  IAZ,  and  IAZ  points  to  the 
inertially  fixed  directions  if  there  are  no  rotations  about  IAX  and  IAy.  Thus,  if,  at  t  =  0  when 
IAX,  IAy<  and  IAZ  are  parallel  to  Xi,  Yj,  and  Zi  axes,  and  if  all  three  control  systems  function 
propertly,  the  platform  would  be  inertially  nonrotating  about  each  of  three  axis  and  the  platform 
becomes  an  inertially  stablized  platform,  even  under  unwanted  disturbances. 

The  actual  three-axis  stable  platform  control  system  is  much  more  complicated.  For  those 
who  are  interested  in  the  further  study  are  referred  to  the  following: 

Inertial  Navigation  Systems 
by:  Charles  Broxmeyer 
McGraw-Hill,  1964 
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Gyroscopic  Theory,  Design  and  Instrumentation 
by:  Wrigley,  Hollister  and  Denhard 
The  MIT  Press,  1969 

Control  System  Design  Analysis  of  Three- Axis  Dynamic  Rate  Table 
by:  Kee  Soon  Chun 

MIT  Dept  of  Aeronautics  and  Astronautics 
(SM  Thesis),  1968 
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16.0  PENDULOUS  INTEGRATING  GYRO  ACCELEROMETER  (PIGA) 


Consider  an  SDOF  RIG  with  a  mass  unbalance  m  located  at  a  distance  r  along  the  S  A  from 
the  center  of  the  gyro  (the  intersection  of  OA  and  SA)  as  shown  in  Figure  16-1. 


FIGURE  16-1.  SDOF  RIG  WITH  A  MASS  UNBALANCE  m 

The  motor-shaft  axis  coincides  with  the  IA  of  the  gu.  The  motor  is  fixed  to  the  base  of  an 
inertially  nonrotating  platform. 

Suppose  the  platform,  such  as  an  IM  of  the  missile,  is  under  acceleration  whose 
nongravitational  portion  (specific  force  fNG)  is  along  the  -X  direction.  This  causes  the  mass  m  to 
experience  the  inertial  reaction  force  toward  +X  direction,  generating  a  torque  about  OA,  that  is, 
about  the  +Y  axis.  If  the  gyro  unit  were  free  to  rotate  about  the  gyro  IA,  the  spin  vector  H  would 
rotate  (precess)  in  the  Y  -  Z  plane  about  the  -X  axis  (Wbc),  trying  to  align  the  H- vector  with  the 
applied  torque  (vector)  about  the  +Y-axis.  The  precession  of  H- vector  is  caused  by  the  inertial 
reaction  force  on  m,  according  to  the  gyroscopic  principle. 
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However,  GU  (whose  input  axis  coincides  with  the  motor  shaft)  is  prevented  from  rotating 
about  its  IA  by  the  fixed  motor  structure.  Therefore,  the  S  A  with  the  mass-unbalance  on  it  has  no 
choice  except  to  move  downwards,  causing  the  SA  to  rotate  about  OA  (+Y  axis)  creating  an  angle 
A  from  the  original  null  position.  This  rotation  is  sensed  by  SG  and  drives  the  motor  shaft 
through  the  servo  motor  control  in  such  a  way  that  the  angle  A  is  nulled  (rebalanced) — in  this  case 

by  rotating  the  motor  shaft  in  the  -X  direction.  When  A  reduces  to  zero,  the  output  of  the  SG 
correspondingly  reduces  to  zero  as  well,  and  the  motor  shaft  rotation  stops. 

For  the  reasons  explained  in  the  previous  sections  on  accelerometers,  fG  (gravitational 
specific  force)  does  not  contribute  to  the  generation  of  torque.  Only  f^c  causes  the  reaction- 
specific  force  fR,  which  generates  the  torque.  The  torque  generated  by  the  inertial  reaction  force 
because  of  f^G  on  the  mass  m  located  at  r  on  the  SA  is  mrfR  =  -mrfNG  with  the  negative  sign 
indicating  that  fR  is  in  the  opposite  direction  from  fNG* 

From  Equation  (16-19)  from  the  Section  13.0  on  the  SDOF  gyro. 


Joa^Aqa  +  CpAoA  -  HSWU  +  mrfNG(l) 


In  (16-1),  the  applied  torque  mrfNG  is  nulled  (rebalanced)  by  the  opposing  gyroscopic 
torque  HsWia  caused  by  the  rotation  of  GU  about  the  negative  input  axis  (about  -X  axis)  by  the 
servo  motor,  thus  causing  the  direction  of  the  torque  generated  by  HsWia  to  point  in  the  -Y 
direction,  which  is  opposite  the  direction  (+Y  direction)  of  the  torque  generated  by  mrfNG* 


Under  steady-state  conditions,  both  PA^  =— Aoa 


and  P2Aoa 


are  zero. 


follows  from  (16-1): 


It 


HsWia  =  -mrfNG 


UW  VIT  jjq 

Writing  ^  ,  fNG=  ^  and  integrating  both  sides  from  0  to  T,  the  computation 

period,  we  have 


That  is,  the  angular  rotation  0  about  IA  accumulated  during  time  T  is  proportional  (scaled 
by  )  to  the  velocity  accumulated  (during  T )  as  a  result  of  the  nongravitational  specific  force 
fNG  applied  in  the  -X  direction. 
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This  device  senses  and  measures  that  part  of  the  velocity,  which  is  the  integration  of 
specific  force  (acceleration)  (caused  by  the  nongravitational  specific  force  f^G)  by  means  of  a  gyro 
with  a  pendulous  mass  deliberately  placed  on  the  SA.  Hence,  it  is  called  a  pendulous  integrating 
gyro  accelerometer  or  PIG  A. 

As  we  have  seen  before  in  the  section  on  SDOF  RIG  (Equation  16-23  of  Section  14.0),  the 
Aoa  of  GE  from  GU  is  related  to  0ia  by 

a  -^e 

AOA  q  dia 

The  Aoa  is  measurable  by  the  signal  generator  attached  to  0A.  Thus,  in  practice,  Vng  may 
be  determined  by  the  summation  of  the  electric  pulses  from  the  signal  generator  attached  on  0A, 
rather  than  by  the  mechanical  summation  of  0ja. 

Unlike  PIPA,  PIGA  does  not  suffer  from  the  voltage  and  current  magnitude  variations. 
Also,  PIGA  is  almost  insensitive  to  the  variations  of  the  magnetic  flux  density  and  permeability. 
This  explains  why  PIGA  replaced  PBPA  in  the  LRBM. 
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17.0  EQUATIONS  OF  MOTION  FOR  INERTIAL  NAVIGATION 


By  Newton’s  second  law: 


eiPjRep  Fg+Fng 


(17-1) 


Since  the  external  forces  have  to  be  either  gravitational  force  (FG)  or  nongravitational  force 
(Fng)  in  a  mutually  exclusive  way,  by  denoting 


***  Q 

f„=—  and  f„„=— — ,  we  may  express  (17-1)  by: 


m 


NG 


m 


PfR™  =- 


I  EP 


m 


1  NG 

m 


=  f„+fi 


NG 


(17-2) 


Equation  (17-2)  is  the  equation  of  motion  to  implement  an  INS  in  the  earth-centered  I- 
ffame.  We  have  previously  found  that  the  onboard  accelerometer  senses  nongravitational  specific 
force  fNG  only,  leaving  out  the  gravitational  specific  force,  fG.  This  necessitates  the  determination 
of  fG  (as  a  function  of  position)  by  the  onboard  computer,  based  on  a  mathematical  model  to 
implement  Newton’s  second  law  by  onboard  instruments  and  a  computer  in  a  self-contained  way 
without  external  help. 

The  block  diagram  for  an  INS  implementation  in  an  I-frame  is  shown  in  Figure  17-1 

below. 
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FIGURE  17-1.  INS  IMPLEMENTATION  IN  AN  I-FRAME 


Now,  to  derive  the  equation  of  motion  in  the  earth-fixed  frame  (that  rotates  with  the  earth 
WRT  the  I-frame  with  angular  velocity  WE),  we  apply  the  theorem  of  Coriolis  to  between  the 
earth-centered  I-frame  and  the  earth-fixed  (and  earth-centered)  E-frame. 

Thus, 


^I^Ep  —  ^E^Ep  Wjg  x  REp 


(17-3) 


It  follows 


pX  =  p.(p,R„) 

=  P,(PEREp  +  WExREp)  using  (17-3) 

=  P[(PePep)  +  P|(Wffi  xR„) .  (17-4) 


Referring  to  the  right-hand  side  of  (17-4),  by  applying  the  theorem  of  Coriolis  to  the  first  term: 
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P.(PeRep)  =  Pe(peRep)+Wie  x(pbrJ 

=  PeRep  +Wffi  x(PeRep)  •  (17-5) 


Noting  that  Wffi  is  constant  and  therefore  PiWje  =  0,  we  have,  from  the  second  term  on  the  right- 
hand  side  of  (17-4) 


P,K  xRj=  (P.W.) xR,tW.x (P,R„) 

=  W„x(P1RB.) 

=  Wjg  x  (PeR-ep  +  Wje  x  Rep) 

=  WE  x  (pEREp)+ WK  x(wffi  xRgp) . 
Substituting  (17-2),  (17-5)  and  (17-6)  into  (17-4): 

■*"  ^ng  =  PeR%  +  2We  X  ^P EREp)  +  Wje  x  x  REp)  . 


(17-6) 


(17-7) 


For  the  case  in  which  the  point  P  is  stationary  on  earth,  REp  is  constant  and,  therefore, 


PeRep  =  °  and  therefore  PE  (PERH>  )=P^REp  =  0. 


For  this  special  case,  (17-7)  reduces  to 


fG-Wffi  x(wffi  x  REp)  =  -fNG 


(17-8) 


As  we  have  previously  discussed  in  the  section  entitled  The  Prototype  Linear 
Accelerometer,  the  accelerometer  output  measures  ■^ng-  Thus,  the  left-hand  side  of  (17-8)  or  fG  - 
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Wffix  (Wffi  x  Rep)  is  the  vector  quantity  (both  direction  and  magnitude),  which  a  plum-bob  on  a 
simple  pendulum  at  rest  on  earth  senses.  We  call  this  vector  quantity  “gravity”  and  denote  it  by  g. 
That  is, 


g£fo- WEx  (wExREp). 


(17-9) 


Note  that  in  (17-9)  g=fG  for  WE  =  0  for  the  nonrotating  earth.  For  the  rotating  earth,  g 
is  not  equal  to  fg  except  at  the  north  and  south  poles,  where  WE  x  REp  =  0. 

Substituting  (17-9)  into  (17-7): 


^ng  +  S  —  PeRep  ■*"  x  (PEREp)  .  (17-10) 

Applying  the  theorem  of  Coriolis  to  (PE  Rep)  between  the  E-frame  and  N-frame  where 
the  N-frame  is  the  local-level  northeast  down  frame, 

P«R«,  =  p.(p,R*) 

=  Pn(PeR^)+WenX(peREp).  (17-11) 


Substituting  (17-11)  into  (17-10),  and  coordinating  (expressing  vector-components)  in  the  N- 
frame: 


~S  "*'1ng  (2We  +  WEN)x(PEREp)  . 


(17-12) 


Equation  (17-12)  is  the  equation  of  motion  for  the  near-earth  (land,  air,  and  sea)  navigation 
implementing  a  local-level  (north,  east,  and  down  (NED)  platform.  It  implies  that  the  velocity 
relative  to  the  earth,  PeRep  and  the  position  relative  to  the  earth,  REp  (which  is  the  integration  of 

PeRep  )  may  be  determined  based  on  the  f"G,  which  is  the  negatives  of  the  outputs  of  three 
accelerometers  along  the  NED  directions  on  the  platform,  which  implement  the  N-frame. 

Referring  to  left-hand  side  (LHS)  of  (17-12),  we  see  that 


[Pn(PbRep)]N=Pn(P.Rep)N 
=  P(PbPep)n 
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Therefore,  since  p  ( • ) = -£•  ( • ) , 

J[p(PERE,)N]dI=Ji(PEREp)',dt 

=  (Mep)N= 

where  Vx,  Vy,  and  Vz  are  respectively  NED  velocities  of  the  vehicle  relative  to  the  earth-fixed 
frame. 

The  vehicle's  position  relative  to  the  earth  (latitude  and  longitude)  is  determined  from  the 
velocity  output  by  integration  in  the  following  way. 


(17-13) 


For  the  spherical  earth,  denoting  latitude  by  <|>  and  longitude  by  X,  and  the  magnitude  of 
Rep  by  R,  we  get  (refer  to  Figure  17-2): 


VN=<j)R  or  <f>  =-gp  so  that 


/  ^dt  =  J<j>dt 


(17-14) 


=  4> 

=  latitude 


and  VE  =  (Rcos<J> )  X  or  X= 


Rcos<|> 


so  that 


f — — — dt=  f  Xdt  =  X  =  longitude. 
J  Rcosd)  J 


Rcos<|) 


(17-15) 
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FIGURE  17-2.  DETERMINATION  OF  LINEAR  VELOCITIES  FROM  THE 
LATITUDE  AND  LONGITUDE  RATES 


Next,  by  denoting  the  magnitude  of  WE  by  Q,  we  have  (refer  to  Figure  17-3): 


( w 

WIEx 

f  £2cos<|>  ^ 

w£  = 

wEy 

= 

0 

IwJ 

£2sin<J>j 

The  negative  sign  in  Wjez  is  necessary  because  the  vertical  component  of  W 
and,  therefore,  opposite  to  the  Z  -  direction,  which  is  downward. 


(17-16) 
is  upward 
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W 

TTffi 


FIGURE  17-3.  DECOMPOSITION  OF  THE  EARTH  RATE  INTO  NED  COMPONENTS 


For  (angular  velocity  of  the  local-level  NED  frame  relative  to  the  earth-fixed  frame), 
the  N-frame  rotates  relative  to  the  E-frame  caused  either  by  latitude  rate  <j>  or  longitude  rate  X. 
Referring  to  Figure  17-4,  the  North  component  is  Xcos<|>  and  the  Z  (down)  component  is  —  Xsincj) , 
because  Xsin<)>  is  in  the  upward  direction,  which  is  opposite  the  positive  Z  (downward)  direction. 
The  north  and  down  components  are  not  influenced  by  <|)  because  <f>  is  perpendicular  to  the 
meridian  plane.  The  east  direction  is  opposite  to  the  <j>  direction  (west).  Therefore,  the  east 
component  -  <f>  is  not  influenced  by  X  because  <j)  is  perpendicular  to  X .  Thus, 


(  Xcos<J>  ^ 

II 

£ 

Weny 

= 

- <i> 

^ — X  sin  <()^ 

(17-17) 
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FIGURE  17-4.  DECOMPOSITION  OF  THE  LONGITUDE  RATE  INTO  NORTH  AND  DOWN  COMPONENTS 

To  maintain  the  platform  (on  which  NED  accelerometers  are  mounted)  in  the  local-level 
(NED)  frame,  it  has  to  be  torqued  by  the  angular  rate  W*  . 

Using  (17-16)  and  (17-17): 

W£  =  W*+W£ 


f  £2cos<|>  ^ 

f  lcos<|> > 

ss 

0 

+ 

-<i> 

£2sin<t>y 

-Xsin<}>^ 
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FIGURE  17-5.  SIGNAL  FLOW  BLOCK  DIAGRAM  OF  LOCAL-LEVEL  N-FRAME 
IMPLEMENTATION  OF  TERRESTRIAL  NAVIGATION 
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FIGURE  17-6.  IMPLEMENTATION  OF  NAVIGATION  EQUATION  (17-12) 


Repeating  Equation  (17-12) : 
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[Pn(PERep)]N=p(p.RbpN)= 

(V  1 

X 

V, 

= 

2  W 

>  > 

V. _ 

y.i 

fO 

fg0 

fN  = 

gN= 

gy 

UJ 

U«; 

w  'l 

w  ^ 

ttENx 

wN  = 

TTffi 

Wffiy 

W^  = 

W 

[wj 

then 

(2W;+W”)x(PsR,p)" 


i  j 

2Wn*+WENx  2Wmy+Wmy 
\  Vy 


2WIEz  +  WENz 


i,  j,  kof  N- frame 


"(2W„  +  WM)V  -  (2W„  +  WHNl)Vy'' 
(2W«,  +  W„)V,  -  (2W„,  +  W„,)V, 
,(2W„  +  W„,,)V,  -(2W„  +  W„,)V,; 


substituting  these  into  (17-12): 

v,  =  f,  +  g,  +  (2W„  +  W„,)V,  -  (2W„  +  W„,)V, 

V,  =  f,  +  g, +  (2^  +  WjV,  -(2W1&  +  W„,)V, 
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V,  =  f.  +  g,+  (2Wmy  +  W„,)V,  -  (2W„  +  W„)V, 
in  which,  referring  to  (17-16)  and  (17-17), 


Wmjt=Qcos<() 

wEy=o 


W1El=-Qsin<j) 


y^Wx=^cos^ 


WENy=-<{) 

W^-lshKj) 
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18.0  DETERMINATION  OF  INITIAL  POSITION  AND  INITIAL  VELOCITY 


In  inertial  navigation,  we  determine  the  acceleration  P*  from 


Pi  Rep(0  ~  "h  ^NG 


(18-1) 


where 


Rep  is  the  displacement  vector  from  the  center  of  the  earth  to  the  point  of  accelerometer 
location,  and  fG  is  obtained  from  the  mathematical  model  as  a  function  of  position  and  fNG  from  the 

accelerometer  measurements.  The  velocity  PjRep  is  obtained  by  integrating  P^Rg,,  that  is, 


PiRepW  —  +  Pi  Rep  (0)  (18-2) 

where  PiRep(O)  is  PiRep  at  t  =  0.  The  position  Rgp  is  obtained  by  integrating  PjRgp  that  is, 

REP(t)  =  JoP1REP(p)dP  +  REP(0)  (18-3) 

where  Rep(0)  is  REP  at  t  =  0. 

Thus,  we  see  that,  to  determine  the  current  position  Rep(0,  we  need  the  initial  velocity 
P!Rep(0)  and  the  initial  position  Rep(0). 


18.1  DETERMINATION  OF  THE  INITIAL  POSITION 
We  start  with 


where 


Rep  — ■  Res  Rsp 


(18-4) 


E  =  Center  of  the  earth 

S  =  Center  of  the  submarine  (ship  or  aircraft)  INS 
P  =  Center  of  the  missile  guidance  system. 

and  the  superscript  I  indicates  that  all  components  are  expressed  in  the  I-frame.  In  (18-4),  the 
submarine's  location  REs  is  either  available  as  an  output  from  its  navigation  system  or  may  be 
computed  from  it. 
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The  position  R*p  of  the  missile  guidance  system,  before  launch,  while  in  the  submarine  or 

aircraft  is  given  in  the  ship-fixed  S-frame.  The  origin  of  the  S-ffame  is  at  the  center  of  the  ship 
navigation  system.  The  x-axis  of  the  frame  is  along  the  longitudinal  axis  on  the  plane  paralleled  to 
the  main  (or  conceptual)  deck  through  the  origin  of  the  frame.  The  y-axis  is  perpendicular,  on  the 
plane,  to  die  x-axis  on  the  plane.  The  z-axis  is  perpendicular  to  both  x  and  y  axes  according  to  the 

right-hand  rule.  To  express  R*p  as  R*p ,  we  have  to  know  C*  so  that, 

1^-sp  =  CSRSP  .  (18-5) 

The  coordinate  transformation  matrix  C*  from  the  S-ffame  to  the  I-frame  is  determined 
from  the  gimbal  angles  of  the  submarine  (or  aircraft)  INS  (SINS)  if  it  implements  the  I-frame.  If  it 
implements  the  navigation  frame  (NED),  C*  is  first  determined  from  the  gimbal  angles  of  SINS, 

and  then  C*  is  obtained  by 


CS=C^.  (18-6) 

The  coordinate  transformation  matrix  from  the  N-frame  to  1-frame  is  given  in  terms  of 
the  submarine's  latitude,  <|>  (an  output  from  the  SINS),  and  the  celestrial  longitude,  Xo  by  (with  the 


derivation  to  follow): 

sin<|>cosXc 

-sinXc 

-cos<t>cosX( 

c 1  - 

-sin<J)sin  Xe 

cosXc 

-cos^sinX, 

^  cos<|> 

0 

— sin  4> 

where 


Xc  =  Xg  +  X 


in  which 

Xg  =  Celestial  longitude  of  Greenwich  meridian 

X  =  Terrestrial  longitude  of  the  submarine  (relative  to  the  Greenwich  meridian) 
Xc  =  Celestial  longitude  of  the  submarine  (relative  to  the  inertially  fixed  reference) 


18.2  DERIVATION  OF  Cj, 

This  derivation  is  shown  here  in  detail  as  a  sample  case  of  determination  of  coordinate 
transformation  matrix  for  the  tutorial  purpose. 


It  turns  out  that  it  is  easier  to  get  C*  first,  and  then  get  C^  by  transposing  C*  or 


18-2 


NSWCDD/MP-95/87 


c^=(cr  )T. 


To  get  (refer  to  Figures  18-la  and  18-lb),  rotate  I-frame  (with  Xi,  Yi,  Z\  axes)  about  Zi  axis 

by  an  angle  of  (celestial  longitude).  Let  us  call  the  resulting  frame  the  A-frame  (with  Xx,  Y  x, 
Zx  axes). 

This  operation  yields  (from  Figures  18-la  and  18-lb) : 


Yx  = 


cos  A.  sin  A.  0YfX,Y 
-sin  A,  cosX,  0  Y, 


o  iJX  z, 


(18-8) 


RA  =  CrR' 


(18-9) 


polar 

axis 


Equatorial 

plane 


FIGURE  18-1  A.  DETERMINATION  OF  C.A 


Xjl  =  XicosXc  +  YisinXc 

X^  =  -XisinX^  +  YjcosXc 
Zx  =  Zj 
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FIGURE  18-1B.  DETERMINATION  OF  C* 


Next,  rotate  A  -  frame  about  the  negative  Y\  axis  by  an  angle  of  <|>  (latitude)  as  shown  in 
Figure  18-2. 

Let  us  call  the  resulting  frame  the  O  -  frame  (with  X$,  and  Z$  axes). 

This  operation  yields: 


N 

4> 

f  cos<()  0  suKfP 

4» 

|XV 

Y* 

= 

0  1  0 

Yx 

lZJ 

^-sin<j>  0  cos<|>. 

A 

(18-10) 


or 


(18-11) 
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Referring  to  Figure  18-3  below,  we  see  that  the  N-ffame  (with  XN  (north),  Yn  (east),  ZN  (down) 

axes)  is  related  to  O  -  frame  (shifted  from  the  center  of  the  earth  to  a  parallel  frame  located  at  the 
vehicle  position  on  the  meridian )  by 


or 


Xn  =  Zq  ;  Yn  =  Y,}, ;  Zn  =  x6 


fx  ) 

N 

O 

O 

N 

yn 

= 

0  1  0 

O 

o 

1 

lZJ 

rN_cnr* 


(18-12) 


(18-13) 


FIGURE  18-3.  DETERMINATION  OF  C* 

$ 

Thus,  from  (18-13),  (18-11),  and  (18-9) 


RN  =  C^R" 


c:c:ra 
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=  C5JC*CiaRi.  -  (18-14) 

Since 

RN  =  CfR1 
we  have 


cfcjcjcf.  (18-15) 

It  follows  using  (18-8),  (18-10),  and  (18-12): 


pH 

o 

o 

N 

(  cos<|>  0  sunjO 

<t> 

( cos  X  sinX  0Y 

CN  = 

'-I 

0  1  0 

0  1  0 

— sinX  cosX  0 

-1  0  0; 

• 

^-sin<|>  0  cos^ 

A 

,  0  0  lj, 

or 


CN  = 
M  — 


r-sin<|>cosXc  -sin<|)sinXc  cos<|> 

-sinXc  cosXc  0 

cos<j>cosX,  — cos<J>sinXc  — sin^ 


(18-16) 


Taking  the  transpose  of  (16)  (exchanging  rows  and  columns)  will  give  Cj,  as  given  in  (18-7). 


18.3  DETERMINATION  OF  INITIAL  VELOCITY  PiRep(O) 

Applying  Coriolis  law  to  Rep  between  the  I-firame  and  the  E-frame  (see  Section  7.0): 


PiRep-  PeRep  +  Wje  x  Rgp.  (18-17) 

Substituting  (18-4)  into  the  right-hand  side  of  (18-17)  for  Rep, 


P iRep  —  Pe(Res  +  Rsp)  +  WE  x  (Res  +  Rsp) 


-PeRes  +  PeRsp  +  We  x  Res  +  Wje  x  Rsp  .  (18-18) 

Referring  to  the  second  term  on  the  right-hand  side  of  (18-18),  applying  the  Coriolis  law  between 
the  E-frame  and  the  S-frame, 
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PeRsp  =  PsRsp  +  Wes  x  Rsp 

=  WES  x  RSP  (18-19) 

because  Rsp  is  fixed  in  the  S-frame,  thus  making  PsRsp  =  0 . 

Substituting  (19  18-)  in  (18-18)  and  noting  that  Wie  +  Wes  =  Wis,  we  have 

PiRep  =  PeRes  +  Wes  x  Rsp  +  Wie  x  Res  +  Wie  x  Rsp 

=  PeRes  +  Wm  x  Res  +  WK  x  RSP .  (1 8-20) 

Evaluating  (18-20)  at  launch  (t  =  0)  and  expressing  it  in  the  I-frame  because  that  is  the  frame  we 
use  for  the  missile  guidance  and  navigation,  we  have 

(P.Rep  f,  =  (PeRe  )',+  (W.  x  Rb  )'  +  (W„  X  Rsp )[  (18-21) 


Equation  (18-16)  shows  that  the  velocity  of  the  missile  at  launch  relative  to  the  center  of  the  earth  in 
the  inertial  space  may  be  determined  in  terms  of  the  values  of  the  three  terms  on  the  right-hand  side 
of  (18-21). 

We  deal  with  them  one  by  one. 

First,  (PeRes  )*  is  the  initial  velocity  of  submarine  or  aircraft  relative  to  the  earth-fixed  frame  with 

components  expressed  in  the  I-frame.  The  term  is  the  initial  value  (at  t  =  0)  of  the 

submarine's  velocity  relative  to  the  earth-fixed  frame  expressed  in  the  N-frame.  It  is  an  output 
from  INS  and  is  transformed  from  the  N-frame  to  the  I-frame  by 

(p.rj:=cuoxp.rj: 

where  C^(0)  is  the  value  of  C*N  att  =  0. 

Next,  for  (WE  x  RB  we  determine  W|  x  R|.  in  E-frame  and  transform  it  to  the  1-frame.  In  the 
E-frame  18- 


(0) 


in  which  |Wie|=W 


where  W  is  the  magnitude  of  the  earth-rate. 
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Thus, 


(W.xRl),- 


1 

0 

|Xft 


J 

0 


k 

w 

Zn 


I E- frame 


(18-24) 


=  -iwy0  +  jwx0 
in  which 


db  _ 

^ES  ~ 


f  V  \ 


vz»y 


is  available  at  an  intermediate  point  of  the  computation  from  INS. 
Then  (WE  x  RB  is  obtained  simply  by 


(w.xrb);=c;(o)(w;xrl)‘  (18.25) 

where  Cg  (0)  is  the  coordinate  transformation  matrix  from  the  E-frame  to  I-frame  at  launch  (at  t  = 

0). 

In  our  case  both  the  I-frame  and  the  E-frame  are  polar-geocentric,  thus  sharing  the  Z  (polar)  axis. 
Thus  Cg  is 


cM= 


( coswx  -sin wx  0^ 
sin  wx  coswx  0 
0  0  1 


(18-26) 


derived  by  rotating  the  E-frame  about  the  polar  axis  relative  to  the  I-frame  by  an  angle  of  Wx, 

where  w  is  the  magnitude  of  the  earth  rate  and  x  is  the  time  elapsed  since  the  two  frames  coincided 
(i.e.,  Xe  axis  with  Xi  axis  and  Ye  axis  with  Ye  axis  while  the  Ze  axis  always  remains  identical 
with  the  Zj  axis). 

Finally  to  evaluate  the  last  term  on  the  R.H.S  of  (18-21),  (WB  x  ,  we  determine  it  first  in 
S-frame  and  then  transform  it  to  the  I-frame  via  the  N-ffame.  That  is, 


(w.xRj;=d(oxs(o)(wjxR;); 


(18-27) 
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In  (18-27),  Rsp  is  simply  the  missile  (its  guidance  system’s)  location  in  the  ship  expressed  in 

S-frame.  is  determined  by  the  guidance  system  by  differentiating  with  respect  to  time  its 

gimbal  angles,  which  are  the  measure  of  the  angular  deviations  between  the  guidance  system's 

inertially  nonrotating  frame  and  the  missile  structure  at  launch  time,  which  is  fixed  to  the  ship.  C ^  ( 0 ) 

is  the  coordinate  transformation  matrix  from  the  S-frame  to  the  N-frame  at  launch,  which  is 
determined  based  on  the  SINs  gimbal  angles  (which  is  the  measure  of  angular  deviation  between 

S-frame  and  N-frame).  C^(0)  is  determined  from  the  latitude  and  celestial  longitude  as  given  in 

(18-7). 
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19.0  CROSS  PRODUCT  STEERING 
(FOR  ROCKET  VEHICLE  GUIDANCE) 


Let  us  denote  the  missile’s  current  velocity  by  VM,  the  required  velocity  (such  as  correlated 
velocity)  by  Vr,  and  the  velocity  to  be  gained  (or  to  go)  by  Vg,  so  that  Vm  becomes  Vr  as  Vq 
goes  to  zero. 

Then, 


Vg  =  Vr  -  Vm  .  (19-1) 

Thus,  if  we  drive  Vq  to  zero  with  thrust,  then  Vm  approaches  Vr.  It  is  intuitively  obvious 
that  Vg  will  decrease  most  effectively  if  thrust  is  directed  parallel  to  the  VG  direction. 

In  fact,  it  turns  out  that,  in  all  current  high-acceleration  rocket  missiles,  if  we  steer  (rotate) 
the  missile’s  thrust  vector  with  the  same  angular  velocity  w  as  the  Vg  vector,  then  Vg  will  be 
driven  to  zero,  thus  guaranteeing  that  the  missile  will  achieve  the  required  velocity  Vr.  To  obtain 
the  mathematical  expression  for  this  w,  denote  the  unit  vector  along  Vq  by  u , 

or 

«=^ 

V„  (19-2) 

By  differenting  (19-3)  WRT  time: 

•  _VoYo-YoV0_v0Yo  YqV0 

X  X  X  (19-4) 


Since  u  is  a  unit  vector,  its  magnitude  remains  to  be  one,  while  its  direction  may  change. 
Thus,  by  denoting  the  angular  velocity  of  the  u-vector  with  respect  to  the  inertial  space  by  w,  we 
have,  as  derived  in  the  section  on  theorem  of  Coriolis, 


u  =  wxu  (19-5) 

where  w  is  the  angular  velocity  vector  of  a  with  its  rotation  axis  perpendicular  to  m  . 


It  follows  from  (19-5): 

u  x  u  =  u  x  (w  x  u)  (19-6) 
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Using  the  vector  identity 

Yj  x(V,  x  Vj)  =  (Y1-Y,)YJ-(V1-VJ)VJ  (19-7) 

we  have  from  (19-6): 

u  x  u  =  (u  •  u)w  -  (u  •  w)u  =  w  (19-8) 

because  u  •  u  =  1  and  u  •  w  =  0  since  w  is  perpendicular  to  u  as  mentioned  above. 


Substituting  (19-3)  and  (19-4)  into  (19-8): 
w  =  uxu 

„Yo./VqYo  YoV 
v0  *{  v*  v02 


(19-9) 


V  V  V  V  V  •  V 

_  -Lo  ^  To-Lo  _  JLo  ^  -Lo  To 

V  V1  V  V2 

To  To  To  To 

The  second  term  of  (19-10)  is  zero  because  Yo  x  Yo  =0.  It  follows: 


(19-10) 


V2  (19-11) 

In  (19-5),  (19-8),  and  (19-11),  we  conclude  that  w  is  proportional  to  the  angular  velocity  of  the 
unit  vector  Vq.  Equation  (19-11)  shows  that  w  may  be  determined  from  Vg  and  VG.  And,  if 

the  rocket’s  thrust  vector  is  steered  with  the  angular  velocity  of  w  as  determined  by  (19-11),  Vg 
will  be  most  effectively  nulled  assuring  Vm  will  become  Vr. 

In  practice,  (19-11)  is  used  in  the  form  of 

w  =  K1(VqxVg)  (19-12) 

where  Ki  is  adjusted  for  good  performance  and  may  or  may  not  be  a  function  of  Vg. 

Denoting  the  missile’s  roll  axis  by  isx,  pitch  axis  by  jay  and  yaw  axis  by  kaz,  in  the 
missile’s  body-fixed  frame,  we  have 
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jB  ^B 

Yo  xVg  =  Vgx  Vgy  Vffl 
Vox  VCT  Vo, 

+J'b(V02V<jx-V<!xVgz)  (19-13) 

+k.(v0xVOT-V0TVax). 

Thus,  the  commanded  angular  velocity  in  pitch  and  yaw  directions  are: 

^ pitch  =  — ^i(YgxYgz  —  YgzYgx)  along  jB  -  axis  (19-14) 

Wy»w  =  Kj  (Vqx Vqy  —  Vqy Vqx  j  along  kB  -  axis  (19-15) 

wroii  does  not  effect  the  direction  of  thrust,  and,  therefore,  is  not  used. 

In  the  LRBM,  w  is  determined  by 

W  =  K1(AtxV0)  (19-16) 

instead  of  (19-12),  where  Ax  is  the  nongravitational  specific  force  ^NG*  which  is  equal  to  the  thrust 
in  a  vacuum.  To  see  that  (19-16)  is  a  modification  of  (19-12),  we  start  with  the  differentiation  of 
(19-1): 

Y,=Y-+Yo  (19-17) 

Since 

Y*  =  g+fNO  =  g+AT  (19-18) 

and  noting  that  VR  is  a  very  slowly  varying  function  and  thus  VR  may  be  approximated  by  zero, 
and  also  noting  that  |Aj.|>>  |  g  |  near  the  deployment  release  point  we  have  from  (19-17)  and 
(19-18):  ~ 

Yg  =  -VM  =  — g  -  AT  =  -At  .  (19-19) 

Substituting  (19-19)  into  (19-12): 

W  =  K(V0x(-AT))  =  K(tXAIxV„)  (19-20) 

which  is  the  same  as  (19-12),  with  K  now  being  a  function  of  time  for  optimum  performance. 
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20.0  DERIVATION  OF  WEIGHTING  MATRIX  W 
FOR  STELLAR  SIGHTING 


20.1  INTRODUCTION 

The  position  and  velocity  of  the  missile  as  determined  by  the  onboard  guidance  system 
during  powered  flight  have  inherent  errors  caused  by  the  errors,  which  exists  in  the  accelerometers 
and  gyros,  as  well  as  those  caused  by  the  launch-time  initial  conditions  associated  with  navigation 
and  fire  control. 

By  means  of  a  stellar  sighting  of  a  predesignated  star  before  to  the  initiation  of  the 
deployment  of  the  reentry  bodies,  the  guidance  system  measures  the  change  in  elevation  AE  and 
change  in  bearing  AB  of  the  star  relative  to  the  guidance  system  platform  axes. 

The  guidance  system  uses  these  measurements  to  make  the  best  estimate,  in  the  least  square 
sense,  of  the  errors  in  position,  velocity,  and  platform  orientation  by  using  a  weighting  matrix 
denoted  as  the  W  matrix.  The  W  matrix  is  predetermined  before  launch,  based  on  the  error 
statistics  of  the  accelerometers  and  gyros  as  well  as  those  of  the  initial  conditions  related  to 
navigation  and  the  fire  control. 


20.2  FORMULATION 


Define  a  state  vector  x,  which  is  to  be  estimated  at  the  time  t.  of  the  stellar  sighting:  bv: 


X  = 


'APl> 

Ap2 

Ap3 

APi 

Ap2 

Ap3 


0! 

e2 

^  03  ; 


(20-1) 


where  Ap,,  Ap2,  and  Ap3  are  the  guidance  system's  position  errors,  Apt  ,  Ap2,  and  Ap3  are 

velocity  errors  and  0i,  82,  and  83  are  platform  tilt  errors,  along  the  X,  Y,  and  Z  axes,  denoted  by 
subscripts  1,  2,  and  3  for  convenience. 

We  denote  a  constant  error  vector  with  some  65  elements  by  e.  Each  element  of  the  e 
vector  is  considered  as  a  random  constant  for  a  given  missile  flight,  but  varies  from  flight-to-flight 
randomly  with  known  standard  deviations  and  mean  values. 
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We  relate  the  state  X  at  the  time  of  sighting  to  e  by  means  of  the  M-matrix,  that  is  X=Me. 
Note  that  X  being  a  (9  by  1)  column  vector  and  e  being  (65  by  1)  column  vector,  the  size  of  the  M- 
matrix  must  be  (9  by  65)  to  be  conformable. 


In  the  stellar  sighting,  the  onboard  instrument  measures  the  elevation  error  AE,  and  the 
bearing  error  AB  (of  the  designated  star)  relative  to  the  guidance  system's  platform  axes.  We 
represent  this  by  the  2x1  measurement  vector  Z  by: 


(20-2) 


We  relate  Z  to  the  state  vector  X  by: 


Z  =  HX  +  v 


(20-3) 


where  v  is  a  2x1  measurement  noise  vector,  and  H  is  called  the  measurement  matrix.  Since  Z  is  a 
2x1  matrix  and  X  is  a  9x1  matrix,  H  must  be  a  2x9  matrix  to  be  conformable. 

Our  objective  is  to  find  the  weighting  matrix  W  such  that  the  estimate  X  of  X  based 
on  the  measurement  Z  is  optimal  (in  a  sense  to  be  defined  later).  That  is,  we  have  to  determine  the 
weighting  matrix  W  such  that  X  in 


X  =  WZ 


(20-4) 


is  optimal. 


It  follows,  substituting  (20-4)  into  (20-5): 

X  =  W(HX  +  v)  (20-5) 

20.3  DETERMINATION  OF  W  MATRIX 

The  state  vector  X  represents  the  true  value  of  the  errors  in  the  guidance  system's 
determination  of  position,  velocity,  and  platform  alignment.  We  seek  the  best  (optimal)  estimate  X 
of  these  guidance  system  errors.  Denoting  by  X  the  error  in  the  process  of  obtaining  the  estimate  X 
of  the  system  error,  we  may  write 

X  =  X  +  X  (20-6) 


It  follows  using  (2)  and  (5)  in  (6)  noting  X  =  Me: 
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X  =  X-X 
=  W(HX  +  v)  -  X 
=WHMe  +  Wv  -  Me 

=  (WH  - 1)  Me  +  Wv  (20-7) 

where  I  is  an  identitiy  matrix.  We  denote  the  covariance  matrix  of  X  by  P*  .  That  is: 

P*AE(XXT)  (20-8) 

Where  E  denotes  the  statistical  expectation  and  the  superscript  T  denotes  a  transpose. 

Substituting  (20-7)  into  (20-8): 

P*  =  e(xXt) 

=  E{[(WH  -  I)Me  +  Wv][(WH  -  I)Me  +  Wvf} 

=  E{[(WH  -  I)Me  +  Wv][eTMT(WH  -  I)T  +  vTWT]} 

=  E{[(WH-I)MeeTMT(WH-I)T]  +  WwTWT} 

with  the  assumption  that 

E(veT)  =  0 
E(evT)  =  0 

This  is  a  reasonable  assumption  because  in  reality,  e  and  v  are  unlikely  to  be  correlated. 
Note  that  the  expectation  E  operates  only  on  the  random  vectors  such  as  v  and  e,  but  not  on  the 
constant  matrices  such  as  W,  M,  or  (WH-I). 

From  (20-10): 

P*  =  E[(WH  -  I)ME(eeT)MT(WH  -  I)T  +  WE(wT)WT] 

=  (WH  -  I)MPeMt(WH  -  I)T  +  WPVWT  (20-11) 

where 


(20-9) 


(20-10) 
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PEAE(eeT)' 

PvAE(wt) 


(20-12) 


We  want  to  determine  W  such  that  the  scalar  sum  of  the  diagonal  elements  of  Px  called  trace  (tr)  of 
square  matrix  Px  is  minimized  because,  generally,  the  diagonal  elements  are  most  significant  and 
predominant.  To  achieve  this,  we  set 

55r«A-°  (2(M3) 


and  we  solve  the  resulting  equation  for  W. 

For  those  who  are  interested  in  a  tutorial  derivation  of  (20-16)  to  (20-22),  see  the  last  paragraph  of 
this  report. 


From  (20-11) 

P*  =  ( WH  -  I)(MPEMTHTWT  -MPeMt)  +  WPv  Wt 

= w(hmpemtht)wt  -  mpemthtwt 

-  WHMPeMt  +  MPeMt  -  WPVWT  (20- 14) 

Now,  since  PE  is  a  covariance  matrix,  PE = PE  and  we  have: 

(hmpemtht)t =(mtht)t(hmpe)t 

=(ht)t(mt)tpet(hmt) 

=  HMPeMtHt  (20-15) 

which  shows  it  is  a  symmetric  matrix. 

According  to  a  matrix  formula ,  denoting  trace  by  tr 

-r^r  tr  (WBWT)= W  (B  +  BT ) 
d  W 

=  2WB  if  B  =  Bt  (20-16) 


Thus: 
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^rtr  W(HMPeMtHt)Wt  =  2W(HMPeMtHt) 
and 


a 

aw 


tr  WPVW  =  2WPV 


since  Pv=Py. 
Now,  by  formula: 


aw 


trAWT  =  A. 


Thus 


^-tr  (MPeMtHt)Wt  =  MPeMtHt 


and  also,  by  formula: 


a 

aw 


trWC  =  CT 


Thus: 


^~trW  (HMPeMt)  =  (HMPeMt)T 

=(mt)t(hmpe)t 

=  M(Pe)t(HM)t 

=  mpemtht 


since  PE  =  Pg. 

Note  that,  since  MPEMT  is  not  a  function  of  W, 


a 

aw 


tr(MPEMT)  =  0. 


(20-17) 

(20-18) 


(20-19) 


(20-20) 


(20-21) 


(20-22) 


(20-23) 
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Thus,  from  (20-14)  using  the  above  equations,  and  setting  the  result  equal  to  zero: 

~~ —  trPjj  =  2W0HMPeMtHt  -  2MPeMtHt  +  2W0PV  =  0  (20-24) 

dW 

with  W  replaced  by  WQ,  the  optimal  value  of  W. 

Solving  (20-24)  for  W0: 

W0(HMPeMtHt  +  Pv)  =  MPeMtHt  (20*25) 

or 

W0  =  MPeMtHt(HMPeMtHt  +  Pv)"!  (20-26) 

or 

W0  =  MPe  (HM)t[((HM)Pe  (HM)t  +  Pv)_1]  (20-27) 

if  the  inverse  exists.  This  is  the  W-matrix  we  are  seeking.  Further  analysis  shows  that  Wo  given 
in  (20-26)  is  indeed  a  minimum  as  required. 

The  optimal  estimate  X  is  obtained  by  (20-5): 

X  =  W0Z 

with  W0  given  by  (20-26) 

For  the  tutorial,  heuristic  derivations  of  the  matrix  formula  used  in  this  section, 
demonstrated  by  means  of  2x2  matrices  for  satisfaction  (not  as  a  rigorous  derivation),  the  readers 
may  refer  to  document  NAVSWC  MP  91-03  "Derivation  of  Discrete  Kalman  Filter  Equations 
using  Heuristic,  Tutorial  Approach"  (Revision  A)  by  Kee  Soon  Chun,  NSWC,  Jan  1991.  Anyone 
interested  in  having  a  copy  may  contact  the  author  at  NSWCDD  (Code  K44). 
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21.0  SCHULER  PERIOD  AND  SCHULER  TUNING 


Consider  a  spherical,  nonrotating  earth.  If  one  stands  still  on  the  earth  and  holds  a  plumb- 
bob,  it  will  point  to  die  center  of  the  earth  in  the  vertical  direction.  If  one  accelerates,  we  know  by 
experience,  the  line  of  the  plumb-bob  will  lag  behind  (from  the  stationary  direction)  because  of  the 
inertial  reaction.  Thus,  we  can  see  that  the  plane  perpendicular  to  the  plumb-bob  line  cannot  serve 
as  the  true  horizontal  plane  because  when  accelerating,  the  plumb-bob  line  deviates  from  the  true 
vertical.  That  is,  the  reference  frame  constructed  based  on  the  plumb-bob  is  not  accurate  because 
the  frame  tilts  under  acceleration. 

Question :  How  can  we  make  the  plumb-bob  follow  the  vertical  even  under  acceleration? 

It  is  intuitively  obvious  that,  if  we  can  make  the  angular  rotation  of  the  plumb-bob  line 
relative  to  its  pivot  point  equal  to  the  angular  rotation  of  the  earth’s  radius  vector  from  the  center  of 
the  earth  to  the  pivot  point,  the  problem  is  solved. 

Now,  consider  a  hypothetical  simple  pendulum  with  the  length  equal  to  the  radius  of  the 
earth  and  the  mass  m  located  at  the  center  of  the  earth.  It  is  obvious  that,  no  matter  how  the  pivot 
accelerates,  the  line  of  the  pendulum  always  points  to  the  center  of  the  earth  along  the  line  of 
vertical. 

From  elementary  physics,  we  know  the  period  T  of  a  simple  pendulum  with  length  L  is: 


T  =  2k 


For  the  pendulum  with  L=R,  the  radius  of  the  earth: 


T  =  2tc  —  =  84  min 

VS 


In  his  classical  paper  (published  in  1923)  on  the  research  of  the  gyrocompass  (North- 
seeking  device),  a  German  scientist,  Max  Schuler  suggested  that  if  we  could  build  a  servo- 
mechanical  control  device  (platform),  which  has  a  natural  period  of  oscillation  of  84  min,  the 
device  will  track  the  earth’s  vertical  under  acceleration  without  oscillation,  if  it  was  initially  aligned 
with  the  earth’s  vertical. 


For  this  reason,  the  84-min  period  is  referred  to  as  the  Schuler  period,  and  the  natural 
frequency  corresponding  to  the  Schuler  period  is  called  the  Schuler  frequency.  That  is, 


Wn  =2tc— =  2ir— 
n  T  2k 


or 
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Wn= 

A  physical  device  (such  as  an  inertial  platform)  that  has  this  arrangement  is  said  to  be 
“Schuler  timed.” 

If  the  platform  is  Schuler  tuned  this  way,  the  angular  velocity  of  the  platform’s  vertical  axis 
is  equal  to  the  angular  velocity  of  the  true  vertical,  if  it  is  intially  aligned.  If  it  is  not  initially 
aligned,  the  platform’s  vertical  axis  will  oscillate  about  the  true  vertical  with  a  natural  frequency  of 
an  84-min.  period  and  with  an  amplitude  of  oscillation  equal  to  the  initial  angle  of  misalignment. 


|-j£-  is  the  Schular  frequency. 
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22.0  CORRELATED  VELOCITY 


22.1  INTRODUCTION 

The  correlated  velocity  is  the  required  velocity  of  a  reentry  vehicle  (RV),  which  will  hit  the 
target  on  the  rotating  earth  by  free-falling  in  a  specified  time  of  flight  (tf).  It  is  the  required  velocity 
correlated  with  the  target  displacement  in  the  inertial  space  caused  by  the  rotation  of  earth  during 
the  specified  TOF  of  the  R V. 

t 

The  specification  of  tf  is  necessary  because  the  target  is  located  on  the  rotating  earth  and 
thus  moves  to  a  new  location  in  the  inertial  space  during  tf,  and  because  we  want  the  RV  to  hit  the 
same  target  at  the  new  location  in  the  inertial  space. 

Another  way  of  looking  at  the  problem  is  how  much  speed  and  in  what  direction  the  RV 
should  have  so  that,  by  free-falling  along  a  elliptic  trajectory,  it  will  intercept  the  target,  which  is 
moving  along  a  circular  orbit  in  the  inertial  space,  while  still  fixed  on  the  earth. 

Referring  to  Figure  22-1,  the  problem  of  finding  the  correlated  velocity  may  be  stated  as 
follows: 

Given:  r(0)  =  ro  at  t=0  of  RV, 

rr(0)  of  target  at  t=0  fixed  on  the  rotating  earth 

tf  time  of  flight  of  free-falling  RV  released  at  t  =  0  until  it  hit  the  target 

Find:  Vc  (Correlated  Velocity  magnitude)  and  yc  (flight  path  angle)  such  that 
r  (tf)  =  rr  at  t  =  tf 

that  is,  such  that  the  position  of  RV  coincide  with  the  position  of  target  on  the 
rotating  earth  at  t  =  tf  or,  RV  intercepts  the  target 

This  problem  is  known  as  Lambert’s  problem  or  the  Lambert  Guidance  problem. 
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Vc ;  Vo 


Z0 


FIGURE  22-1.  CORRELATED  VELOCITY 

The  derivation  of  equation  to  determine  the  correlated  velocity,  in  the  approach  adopted  in 
this  report,  in  terms  of  ro  and  ir(tf)  is  based  on  two  key  equations. 

The  first  equation  is  the  expression  for  the  conservation  of  angular  momentum  along  the 
line  normal  to  the  plane  containing  the  RV’s  trajectory.  The  second  equation  is  the  equation  of 
motion  for  the  RV  derivable  from  Newton’s  second  law  of  motion  along  the  radial  direction  of  the 
RV  from  the  center  of  the  earth. 
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where  Pj(-)  represents  )  in  inertial  space  and  M  is  the  mass  of  the  earth  that  is  assumed  to 
be  concentrated  at  its  center.  Using  Ptr=£,  P^r=r  and  |i  =  GM  in  (22-1): 


multiplying  both  sides  of  (22-2)  by  r  x  (vector  cross-product): 

u. 

rxr  = — -rxr 

-  -  3  "  " 


(22-2) 


(22-3) 


Noting  that  rxr  =  0,  we  have  from  (22-3): 

I  xr=0 


Since,  by  simple  differentiation 


d  .  .  _  .. 

-rrxr=rxr+rxr 
dt  ------ 


and  since  £  x  £=Q,  we  have  using  (22-4): 


d_ 

dt 


(rx£)  =  0 


or 


r  xf=h= constant  vector. 


(22-4) 


(22-5) 


(22-6) 


(22-7) 


This  means  for  anytime  t, 

l(0)x£(0)=i(t)x£(t)  (22-8) 

=  h 

To  proceed  further,  we  define  a  new  rotating  frame  of  reference  that  we  call  the  T-frame. 
Referring  to  Figure  22-1,  the  x-axis  is  along  the  r(t)  direction  and  the  y-axis  is  along  the  line 
perpendicular  (in  counter-clockwise  direction)  to  r(t)  on  the  trajectory  plane.  The  z-axis  is  normal 
to  the  x-y  plane  (trajectory  plane)  with  the  origin  at  the  center  of  the  earth  sharing  the  origin  with 
the  earth-centered  I-frame.  The  T-frame  rotates  about  the  z-axis  with  the  angular  velocity  Wir=S. 
of  the  T-frame  relative  to  the  earth  centered  1-frame. 

Since  the  h  vector  is  by  definition  of  cross  product,  perpendicular  to  both  r  and  £ ,  it  is 
perpendicular  to  the  plane  formed  by  iand£.  Since  h  is  constant,  the  plane  containing  i  and  £ 
is  constant  as  well  for  all  t,  which  means  that  the  trajectory  of  RV  remain  on  a  fixed  plane. 
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Applying  the  law  of  Coriolis  to  r(t)  between  the  I-frame  and  the  T-ffame  (refer  to  the 
section  on  the  law  of  Coriolis  for  derivation): 

PIr(t)  =  PTr(t)+Wn.xr(t)  (22-9) 


where 


I 

Pi  (  )  ”  dt  ( ' )  —  ^me  rate  change  of  vector  (•)  as  observed  in  the  I-frame. 


dT 

PT  (  ) = "dt"  (  ^ —  hme  rate  of  change  of  vector  (•)  as  observed  in  the  T-ffame. 


In  the  T-frame,  according  to  our  definition:  (Note:  superscript  T  should  not  be  confused  with 
transpose.) 


(22-10) 


(22-11) 


so  that 


i  j  k  /0  )T 

Wn.xr(t)=  0  0  0  =  r0  expressed  in  T-frame . 
r  0  0  x  k  0  > 


(22-12) 
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Substituting  (22-10)  and  (22-12)  into  (22-9) : 


T 

0  ' 

T 

'r  \ 

0 

+ 

r0 

= 

r0 

,0, 

i  0  ; 

,0, 

with  components  expressed  in  the  trajectory  frame  (T-frame). 
It  follows,  expressed  in  the  T-frame: 


i  j  k 

r  0 

hT=rTxfT= 

r  0  0 

Z= 

0 

rr0  0 

T 

/*> 

(22-13) 


(22-14) 


According  to  (22-14),  the  angular  momentum  viewed  from  the  trajectory  plane  is  pointed 
along  the  T-frame’s  z-axis  as  it  should  be  and  is  equal  to 


h=r20=r 


2d0 
dt 


From  (22-15),  we  get 


(22-15) 


(22-16) 


(22-16)  is  the  first  of  two  key  equations  necessary  to  derive  the  equation  to  determine  the  correlated 
velocity  in  terms  of  known  quantities. 

EQUATION  OF  MOTION — along  the  x-axis  (radial  direction)  of  the  trajectory  (T)  frame. 
Operating  on  both  sides  of  (22-9)  by  Pi : 


PI(PIr)=PI(PTr)+PI(Wirxr).  (22-17) 

Referring  to  the  right-hand  side  of  (22-17),  treating  Pxr  and  W^xr  as  vectors  and  applying 
Coriolis  law: 
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Pi(V)=PT(PTO  +  Wirx(PTr) 

=  P"r+Wirx(PTr) 


and 


P,(W„xr)=PT(Wirxr)+WItx(WIrxr) 

=  (PIWff)xr+Wlrx(PTr)+WIrx(Wlrxr). 

Substituting  (22-18)  and  (22-19)  into  (22-17): 

P?  r-  P^+  (PT W„)  x  r + 2W„  x  (P^) + W„  x  (Wff  x  r) . 


Noting  that  0  is  along  the  z-axis  and  thus  Wn.= 


O' 

T 

(x 

0 

,  andr= 

0 

0, 

,0; 

,  we  have 


P*r=PT(PTr)=P1 


'i' 

0 

loj 


\T 

0 

loj 


P  W  =P 

X'fTTjj.  *  T 


O' 

f°f 

0 

— 

0 

fit 

(PTWff)xr= 


i  j  k 

0  0  0 
r  0  0 


'0  V 

r0 

loj 


Wirx(PTr)= 


i  j  k 

f0 

0  0  0 

- 

r0 

too 

T 

,0, 

and  using  (22-12) 


(22-18) 


(22-19) 


(22-20) 


(22-21) 


(22-22) 


(22-23) 


(22-24) 


22-6 


NSWCDD/MP-95/87 


W,r*(Wir*r)= 


i  j  k 

-r62' 

0  0  9 

- 

0 

0  r6  0 

T 

,  0  ; 

(22-25) 


Substituting  (22-21),  (22-23),  (22-24),  and  (22-25)  into  (22-20): 


'f ' 

0  ' 

m 

-re2^ 

0 

+ 

r9 

+ 

•o 

+ 

0 

,0; 

1 0  ; 

lo  J 

,  0  , 

(22-26) 


with  components  expressed  in  the  T-frame.  It  follows  that  the  x-component  (in  the  radial  direction 
of  (22-26)  is: 


(p?r)> 

and  the  y-component  is 

(Pjr)^=r0  +  2re 

Since  the  gravitational  force  is  along  the  radial  or  x-axis: 


±LrT_iLJL_ 

r3  -  r>  r 


7  it' 
r2 
0 

l  0  ) 


r  . 


where  ^  is  a  unit  vector  along  r  direction. 


From  (22-2),  (22-27),  and  (22-29): 


(22-27) 


(22-28) 


(22-29) 


(22-30) 
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in  the  radial  direction  or  along  the  x-axis. 

And,  from  (22-2),  (22-28),  and  (22-29): 

r0  +  2r  0  =  0  (22-31) 

in  the  direction  perpendicular  to  the  radial  direction,  or  along  the  y-axis  direction. 

Equation  (22-30)  is  the  equation  of  motion  of  the  RV  along  the  radial  direction,  the  second 
key  equation  we  are  looking  for. 

Referring  to  (22-31): 

r@  +  2rG = —  (r2§  +  2rr0 ) 

=— -j-  (r20  )=0.  (22-32) 

r  dt  '  ' 


From  (22-32): 

Integrating,  we  have  r20  =  constant.  Denoting  this  constant  by  h,  we  have 

r20  =  h 


(22-33) 

(22-34) 


which  is  identical  with  (22- 15), thus  reaffirming  the  conservation  of  angular  momentum. 


Defining  a  new  variable  u  by 

and  substituting  in  (22-34): 
Now 


1 

U4-  or  r= 
r 


\ 

u 


0  =  hu2 


dr  _  dr  d0 
dt  -  d0  dt 


9  d0  9  d0 


(22-35) 

(22-36) 


Using  (22-36)  in  the  above  equation  : 


Also, 


dr  dr  d0  ^  dr 
dt  "  d0  dt  -W  d0 


(22-37) 
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Using  (22-36)  and  (22-37)  in  the  above  equation: 


Finally,  using  (22-35),  (22-36),  and  (22-38)  in  (22-30): 
which  is  a  linear  differential  equation  equivalent  to  (22-30). 


(22-38) 


(22-39) 


22.3  CORRELATED  VELOCITY 

This  section  is  a  detailed  version  of  the  method  of  determining  the  correlated  velocity 
described  in  Reference  (7),with  the  gaps  and  details  filled  in  by  the  author  of  this  report 

There  are  several  approaches  to  the  solution  of  the  problem.  Although  the  method  adopted 
here  is  not  efficient  in  computation  time,  it  was  chosen  because  it  is,  from  the  author’s  viewpoint, 
the  least  difficult  and  least  prone  to  losing  the  insights  into  what  is  going  on.  The  most  efficient 
method  may  be  somewhat  difficult  in  providing  insight  for  beginners.  Those  readers,  who  are 
interested  in  further  study  of  the  subject,  are  referred  to  an  excellent  book  by  Richard  H.  Battin 
entitled  “  An  Introduction  to  the  Mathematics  and  Methods  of  Astrodynamics”  (AIAA  Education 
Series,  1987.),  which  must  be  counted  among  the  great  classics  in  scientific  literature  as  the  Editor- 
in-Chief  of  the  series  noted. 

Recalling  that  (Figure  22-1) 


h=I0x£=r0V0sinY 


(22-40) 


and  referring  to  the  —  term  in  (22-39),  we  may  write 
h 


h2  r2V2sin2y 


r0sin  y 


1 

Ar0sin2Y 


(22-41) 
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where 


Substituting  (22-41)  into  (22-39) : 


d2u(e) 

d0 


+u(e)  J  . 

Ar0sm  y 


(22-42) 


(22-43) 


Referring  to  Figure  22- 1 ,  assuming  6  =  0  when  t  =  0  and  r  ( 0 )  =  r0  ,  the  initial  conditions 
necessary  for  the  solution  of  (22-43)  are: 


U(°)=7M  (22_44) 

and  referring  to  Figure  22-2, 


d  /A.  d  1 
-5eu(e)=de7(9) 
__  1  dr 

T2  d0 

1  dr 
r  rd0 
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= — -coty 


Thus,  since  r  =  ro  when  0  =  0q  : 


du(0) 

d0 


_ _ 1  cosy 

e=o~  ro  siny 


(22-45) 


The  complete  solution  of  (22-43)  is  (as  shown  in  the  appendix  at  the  end  of  this  section): 

1  1  ( l-cos0  sin(y-0) 

I  .  ,  • _ 


u(0)  = 


r(9)  r0 


{  Xsin2y  siny  j 


(22-46) 


The  validy  of  (22-46)  as  the  solution  of  (22-43)  may  be  confirmed  by  differentiating  u  (0 )  W.R.T 

e, 


and 


du (0)  1 

sin0 

cos(y-0) ) 

d©  "r„ 

^sin2y 

siny  I 

d2u(0)  1 

COS0 

sin(y-0) ' 

1  ° 
l-t 

II 

ts 

CD 

k  Xsin2y 

siny  . 

(22-47) 


(22-48) 


and  substituting  (22-46)  and  (22-48)  into  (22-43). 

When  0  reaches  <J>  (range  angle),  r  becomes  rT,  the  target  vector. 
Evaluating  (22-46)  at  this  boundry  condition. 


'  1  -  cost))  ^  sin  (y  - 
X.sin2y  siny 


(22-49) 


Solving  (22-49)  for  X  and  recalling  the  definition  of  (22-42) 


X  = 


1  -  COS([) 


(h\ 

\Tr) 


sin2y  +  sin  ( <|>  -  y )  siny 


(22-50) 


Solving  (50)  for  Vq: 
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\{\i) 

1  -  COS(|)  ] 

lro  J 

fr°l 
lrT  J 

sin2y  +  sin  (<(>  -y )  siny 

(22-51) 


There  are  many  combinations  of  Vo  and  y  (magnitude  and  direction  of  correlated  velocity) 
which  satisfy  (22-51)  for  given  values  of  ro,  ry  and  <|> .  To  determine  the  unique  values  of  Vo  and 
y ,  we  have  to  specify  the  time  of  (free  fall)  flight,  which  is  shown  next 

Referring  to  (22-16),  integrating  dt  from  0  to  tf  (time  of  flight)  and  d0  from  0  to  <{>  (range 
angle),  we  have 


<t> 

JV(9)de 

0 


(22-52) 


Substituting  for  r(0)  from  (22-46),  and  for  h  from  (22-40): 


r  *r 
t  0  f 

'  1-COS0  ' 

+  j'sin(y-0)  \ 

-2 

d0 

f  V0sinyJ 

,  ^Sin2y  , 

{  siny  j 

(22-53) 


with  X  given  by  (22-50),  and  Vo  given  by  (22-51). 

Equation  (22-53)  shows  that  for  the  given  time  of  flight  tf  and  range  angle  <(),  the  Vo  and  y, 
which  satisfy  (22-53)  becomes  the  correlated  velocity  Vc  =  Vo  with  the  flight  pass  angle  yc  =  y  at 
t  =  0  or  at  deployment. 

The  integrand  of  (22-53)  may  be  written 


l-cos0  sin(y-0) 

-2 

1  cos0  siny  cos0-  cosy  sin0 

-2 

Xsin2y  siny 

Xsin2y  Xsin2y  s^nT 

1 

,  + 

f1-  \  1 

cos0-cotysin0 

Xsin2y 

Xsm  y  J 

' 
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=  [a+  bcosG  +  csinG  ]~2 


(22-54) 


where 


1 

Xsin2y 
b=  1  -  a 
c=-coty 


Using  (22-54)  in  (22-53) : 


V0siny 


♦ 

f  d8 

q  (a+bcos0  + csinG)2 


(22-55) 


The  integration  may  be  performed,  resulting  in  a  closed  form  solution  based  on  the  formula  table. 
After  some  algebra,  the  final  algebraic  equation  for  tf  may  be  written : 


t_  ro 

r 

( 1  -  cos<|> )  coty  +  ( 1  -  X )  sin<|) 

2siny 

((2/X)- 

1)*  1' 

f  V0siny  • 

(2-X) 

(r  ) 

Ao 

T 

X 

(f-.l 

tan 

3/2 

siny  cot 

[♦1 

U  J 

> 

-cosy 

(22-56) 


Referring  to  (22-56)  and  equivalently  to  (22-53),  we  note  that  tf  becomes  infinite  at  y  =  0, 

since  siny  is  in  the  denominator.  However,  for  the  current  LRBM  applications,  y  never  becomes 
zero  as  explained  below. 

Multiplying  both  sides  of  (22-49)  by  Xsin2y ,  we  have 

—  Xsin2y  =  1  -  cos<|>  +  Xsiny  sin  (y  -  <|> )  (22-57) 

rT 

It  follows  from  setting  y=  0  in  (22-57)  that  we  will  have  cos<(>  =1  or  the  range  angle  <b  = 
0.  This  means  that  the  current  position  vector  is  r^  colinear  with  the  target  vector  rx  which  never 
happens  in  the  existing  LRBM  practices.  Of  course,  this  is  intuitively  obvious  without 
mathematics. 

The  objective  of  the  following  is  to  verify  Equation  (51)  for  the  special  case  of  a  circular 
orbit  (satellite  injection  into  a  circular  orbit). 
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For  a  special  case  of  a  circular  orbit,  1 I  =  I  rx  I  and  7=90°.  From  the  range  angle  <|> 
introductory  physics,  we  know,  for  a  circular  orbit  of  r0  that 


It  follows: 


(22-58) 


(22-59) 


Evaluating  (22-51)  with  r0  =  rT  and  y— 90*  and  noting  that 


sin  (<{>  -  y  )  =  sin<|>  cos90  °  -  cos<|)  sin90  ° 
=-cos<|) 


(22-60) 


we  have 


p.  1  -  cos<ft 
r0  1  —  cos<{> 


H 


(22-61) 


Thus,  verifying  the  validity  of  (22-51)  for  the  special  case  of  a  circular  orbit. 
Now,  for  a  circular  orbit  of  period  TP,  V0T  =  27tr0.  That  is 

_  2jct0 

T  = - 

p  vo 

For  a  range  angle  <|>  and  the  corresponding  time  of  flight  tf  in  a  circular  orbit : 


Using  (22-62)  in  (22-63), 


Jl=!_ 

Tp  2ti 


(22-62) 


(22-63) 
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<K 


tf==  V« 


Also  for  a  circular  orbit,  using  (22-59)  in  (22-42) : 


Xa 


r  V2 
*o  vo 


(22-64) 


=  1 


Using  7  =  90%  X=1  in  (22-53)  and  noting 

sin  ( y  -  0 )  =  sin90  °  cos0  -  cos90  °  sin0 
=  COS0 

we  have 


T  =—  f 

f  V0  J 


1-COS0  +  COS0 


d0 


(22-65) 


(22-66) 


vo 

which  is  equal  to  (64). 

Thus  verifying  the  validity  of  (53)  for  the  special  case  of  a  circular  orbit. 
22.4  TRAJECTORY:  A  SEGMENT  OF  ELLIPTIC  ORBIT 


(22-67) 


The  general  solution  of  (22-39)  is  given  by 


u(0)=7^y=^-(1+ecose) 


This  may  be  verified  by  differentiating  (22-68) : 

du(0) 

d0 


u, 

—  ,  esin0 
h2 


(22-68) 


(22-69) 
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d2u(6) 

de2 


u. 

— -ecos0 
h2 


(22-70) 


and  substituting  (22-68)  and  (22-70)  into  (22-39),  and  obtaining  the  right-hand  side  of  (22-39). 
Rewriting  (22-68)  in  a  more  conventional  form: 


h2/u 
1  +  ecosO 


(22-71) 


In  (22-71),  —  =  a(l-e2)  as  shown  in  Appendix  B.  Thus,  (22-71)  becomes  : 

r 


r(8)  = 


a(i-ez) 
1  +  ecos0 


(22-72) 


which  is  the  equation  for  the  conic  sections  (see  for  instance,  “Calculus  and  Analytic  Geometry”  by 
Thomas,  Addison  Wesley)  where  e  represents  the  eccentricity  of  the  conic  sections,  and  a  is  the 
semi-major  axis  in  the  case  of  ellipses.  It  is  well  known  that  the  conic  section  becomes  a  circle  if  e 
=  0,  an  ellipse  if  (keel,  a  parabola  if  e=l  and  a  hyperbola  if  e>l. 

Now,  solving  (46)  for  r(0): 


rAsin2y 

r(0)=l  +  [>.sinY(sinY-0)-cos0]  (22'73) 

which  is  in  the  form  of  (22-72). 

From  (22-72)  and  (22-73),  we  can  see  that  the  relationships  between  e  and  X  may  be  determined. 
It  turns  out,  after  some  algebra,  that : 

71 

X  =  1  corresponds  to  e  =  0  for  circle  ( with  y-—). 

0  <  X  <  2  corresponds  to  0  <  e  <  1  for  ellipse 

X  =  2  corresponds  to  e  =  1  for  parabola 

X  >  2  corresponds  to  e  >  1  for  hyperbola 

In  applications  for  LRBM,  the  combinations  of  ro,  vo  and  y  is  such  that  X  is  always  in  the  range 

of  0  <  X  <  2,  which  corresponds  to  0  <  e  <  1,  thus  assuring  the  free  falling  ballistic  trajectory  is 
always  a  segment  of  an  elliptic  orbit 
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22.5  COMPUTATION 

The  key  equation  for  the  computation  of  the  correlated  velocity  is  the  Equation  (22-56). 
What  we  want  is  the  magnitude  of  the  velocity  V  o  and  its  direction  y  which  satisfy  (22-56)  for  a 
specified  time  of  flight  tf,  with  the  given  values  of  the  initial  release  or  deployment  position  ro,  the 

target  position  rr,  and  the  range  angle  <|>  between  £,  and  rr .  Since  by  the  defintion  of  the  dot 
product 

r0  •  rT=r0rTcos<|) 

we  can  determine  <|>  from 


<(>  =  cos 


■  i  lo  ’ 

r0  rT 


Note  that  X  in 


(22-56)  is  given  by  the  definition  : 


M- 


(22-42) 


This  implies  that  Vo  appears  in  several  places  in  (22-56)  either  explicity  or  implicity. 

For  the  determination  of  Vo  and  y  by  a  numerical  method,  we  start  with  an  initial  trial  guess  of  y, 
say  45*,  and  get  the  value  of  Vo  corresponding  to  this  value  of  y  by  means  of  equation  (22-51). 

Next,  substitute  the  values  of  y  and  Vo  obtained  above  to  see  if  the  right-hand  side  of  (22-56) 
matches  the  LHS  of  (22-56),  which  is  the  specified  value  of  the  time-of-flight  tf.  Otherwise,  we 

continue  the  process  by  adjusting  the  trial  value  of  y,  which  yields  a  new  value  of  V  o,  until  the 
difference  between  the  LHS  and  the  RHS  of  (22-56)  falls  within  the  acceptable  limit.  The  Vo  and 

y,  which  satisfy  this  condition  is  the  correlated  velocity  we  are  looking  for  under  the  specified 
conditions. 
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APPENDIX  A 

SOLUTION  OF  EQUATION  (22-43) 
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Repeating  (43): 


d2u(e) 

de 


+u(e)= 


1 

XV0sin2y 


Thus,  with  — 
do 


=  D» 


(d2+  1  )u= - — —  =  K 

Xr0sin2y 

Let  Uh  denote  the  homogeneous  solution  and  Up  particular  solution  of  (A-l). 
homogeneous  solution: 


(d2+i)=o 

(D+i)  (D-i)  =  0 


D=i  D=— i 


Thus, 

For  the  particular  solution: 

Assume 

Up  =  A 

then 


uh  =  CjCOS0  +  C2  sin0 


Substituting  into  (A-l) : 

0  + A  =  K 

so,  the  complete  solution  u = uh  +  up  is 
u=  CjCosO+CjSinO+k 
Repeating  (22-44) 


A-3 


(A-l) 

Then  for  the 
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at  0=0 


Substituting  (A-2)  evaluated  at  0  =  0  into  (A-3): 
— =C,+k 


or 


(A-3) 


Differentiating  (A-2): 
du 

— =-Cj  sin0  +  C2cos0 . 

Repeating  (22-45): 

_  1  cosy 

6=0  ~  ro  siny 


du 

d0 


Substituting  (A-6)  into  (A-5)  and  evaluating  at  0=0: 

^  _  1  cosy 
2  r0  siny 


(A-4) 


(A-5) 


(A-6) 


(A-7) 


Substituting  (A-4)  and  (A-7)  into  (A-2)  with  K= - — 

ta0sin  y 


u= 


r_l_ 

1 

ro 

ta0sin2y 

.  1  cosy  .  1 

cos0 - : —  sm0+- 


smy 


Xr0sin2y 


1 

Xsin2y 


cosy 

siny 


sin0  h - — 

Xsiny 
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<x> 

</> 

g 

1 

rH 

|  sinycosG-cosysinG 

ll  > 

siny 

r  1-cosQ  ^  sin(y-G) ' 
Xsin2y  siny 


wMch  is  identical  with  (22-46)  in  the  main  text 
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APPENDIX  B 
EQUATION  OF  ELLIPSE 
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FIGURE  B-l.  ELLIPSE 

The  equation  of  an  ellipse  in  rectangular  coordinate  is,  as  we  all  know,  from  analytic  geometry: 


x 


2 


(B-l) 


The  equivalent  form  of  an  ellipse  in  polar  coordinates  with  r(0)  measured  from  one  of  the  focal 
points  (f)  is  given  by : 


r(0)  = 


a(l~e2) 
1  +  ecosQ 


(B-2) 


See,  for  instance,  Thomas:  Calculus  and  Analytic  Geometry)  where  the  eccentricity  e  is  the  ratio 
of  distance  (of)  to  (oa)  (i.e.,  e  =  -^7  }  The  e  has  to  be  positive  and  less  than  1  (i.e.,  0  <  e  <  1)  to 


oa 


have  geometrical  shape  of  ellipse. 


The  semimajor  axis  (a)  may  be  chosen  arbitrarily  as  long  as  it  is  a  positive  real  number.  So,  we 
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pick  a  number  such  that 


(B-3) 


where 


,  rxmr 

h= - =rxr 

m 


is  the  angular  momentum  per  unit  mass  and 


p=GME 

in  which  G  is  the  universal  gravitational  constant,  and  Mg  is  the  mass  of  the  earth. 

Since  0  <  e  <  1  for  an  ellipse,  e2  <  1.  Thus,  the  RHS  of  equation  (B-3)  assures  a  >  0  because  h2  > 
0  and  |i  >  0. 

Substituting  (B-3)  into  (B-2),  we  have 


r(9)  = 


h2/n 

1  +  ecosG 


which  is  identical  with  (22-71). 
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